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A UNIFIED SPECIAL PERTURBATION MODEL FOR THE 
MOTION OF THE EARTH-MOON SYSTEM 

By 

W.J. Breedlove, Jr. 1 
SUMMARY 


This report contains a theoretical development ' of the equations 
of motion governing the Earth-Moon system. The Earth and Moon are 
treated as finite rigid bodies and a mutual potential is utilized. 
The Sun and remaining planets are treated as particles. Relativ- 
istic, non-rigid, and dissipative effects are not included. 

The translational and rotational motion of the Earth and 
Moon are derived in a fully coupled set of equations. Euler 
parameters are used to model the rotational motions. 

The mathematical model developed herein is intended for use 
with data analysis software to estimate physical parameters of 
the Earth -Moon system using primarily LURE type data. 

The Appendix contains two program listings. Program ANEAM0 
computes the translational/rotational motion of the Earth and 
Moon from analytical solutions. Program RIGEM numerically inte- 
grates the fully coupled motions as described above. 


l 


Associate Professor of Engineering, School of Engineering, Old 
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INTRODUCTION 


The Lunar Laser Ranging Experiment (LURE) , (ref. 1) has 
resulted in the placement of three ranging retroref lectors on 
the Moon. A series of measurements of the distance of these 
retroref lectors from several Earth-based observatories began in 
August 1969. These measurements, at present, allow the deter- 
mination of the distance to the Moon with an accuracy of ± 8 cm. 

A resolution of ± 2 to 3 cm is expected within the next few years. 
Overall, a ± 10-cm accuracy over a 10-year period will soon be 
available. 

The LURE data, in combination with other data types, can be 
used to determine parameters related to the internal composition 
of the Earth and Moon (ref. 24) . Also, checks of current gravi- 
tational theories may be made (ref. 1) . For example, data 
accuracies of ± 3 cm would make feasible the determination of 
the following parameters (refs. 1 and 2): 

A. Geometrical and Orbital Parameters ; 

station coordinates 
retroref lector coordinates 

Earth and Moon orbital constants of integration 

B. Geophysical Parameters; 

station drift 
polar wobble 
rotation variations 
Earth tide 

universal time determination 
orbital acceleration 

C. Selenophysical Parameters; 

physical librations 
free librations 
Moon tide 
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D. Systematic Error Sources; 
fixed bias 

zenith-distance bias 
arbitrary periodic biases. . 

The original mathematical model (October 1973) of the LURE 
team (ref. 1) involved (l) the numerical integration of the Moon 
and major planets as point masses including relativistic effects, 

(2) the utilization of an analytical lunar physical libra tion 
theory based on Eckhardt's work, plus certain additive and 
planetary terms from the Improved Lunar Ephemeris and 3rd and 
4th order terms in the lunar potential (refs. 19 and 23), -and 

(3) determining the angular position and pole of the Earth from 
BIH data. Earth tides and dissipative effects were not modeled 
(ref. 2) . Lunar orbital-rotational coupling was not fully 
modeled (ref. 22). Finally, the BIH data imposes limits on the 
accuracy achieveable from this model. 

The model described above reached its current state by 
appending additional effects to existing models. For example, 
Eckhardt's original libration theory did not include the 3rd and 
4th order terms in the lunar potential. Also, the additive. and 
planetary terms were appended to this original theory. 

This model provided (October 1973) rms residuals in range 
of ± 3 meters. An improved libration theory would considerably 
reduce this value. The LURE team suggested that a numerical 
integration of Euler's equations holds promise for future gains 
in accuracy for the rotational motion of the Moon (ref. 19). 

The above residuals imply the existence of unmodeled effects 
or modeling inaccuracies (ref. 2) . 

The determination of geometrical and orbital parameters, 
geophysical parameters, selenophysical parameters, and systematic 
error sources to an accuracy compatible with the observational 
accuracy thus awaits the development of a rigorous model of the 
Earth-Moon system. Previous models have been attempted in a 
piecewise fashion. Thus, there is a need for a consistent 
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theoretical, mathematical model incorporating Earth and Moon 
rotational, translational, and deformational motions in a 
coupled sense. This model should allow for an inhomogeneous 
Earth and Moon, dissipation effect, general relativity effects, 
and planetary perturbations. Secondly, there is a need for a 
"special" numerical model incorporating pertinent effects from 
the above theoretical model to be used in the parameter estimation 
process. 

The "special" model envisioned at this point (although subse- 
quent investigations may modify it) is of the following form. 

1. Treat Sun and planets as, perturbing point masses, 

2. treat Moon as a tidally deformed body, 

3. treat Earth as a tidally deformed body, 

4. consider the coupled orbital-rotational motions of the 
Earth and Moon (ref. 22) , and 

5. consider the effects of relativity (refs. 25 and 26). 

The governing equations of motion for this "special" model 
are to be numerically integrated (refs. 4 and 9). 

The appropriate numerical integration routine to be used in - 
this model should be investigated. One candidate is an extremely 
accurate Cowell type routine used by Oesterwinter and Cohen in a 
determination of planetary masses (ref. 7) . This routine has 
been developed over a period of years by Cohen and Hubbard and 
has been used primarily for solutions to the planetary n-body 
problem. This scheme is based on the use of a 16th-order set 
of predictor-corrector formulae for integrating accelerations. 
Herrick (ref. 27) also espouses the use of numerical integration 
schemes that integrate accelerations directly. 

GENERAL SYMBOLS AND NOMENCLATURE 
vector 

universal gravitational constant 
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r i 



u 
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first and second time derivatives 

radius vector from solar system barycenter to mass 
center of body i (i = 1,11) 

radius vector from Sun to mass center of body i 
(i = 1,11) ■ 

gradient operator with respect to coordinates of 
mass center of body j (j = 1,11) 

work function 

mass of body j {j = 1,11) 

Euler parameters locating Earth body axes with 
respect to Earth reference axes (i = 0,1, 2, 3) 

absolute angular velocity components of Earth 
resolved along body axes (i = 1,2,3) 

absolute angular velocity components of Moon 
resolved along lunar body axes (i = 1,2,3) 

Euler parameters locating lunar body axes with 
respect to lunar reference axes (i = 0,1, 2,3) 

components of spherical polar coordinate system 

vector-column 

matrix 

transposed vector -row 
translating reference frame 
vector-row 
dyadic 


PHYSICAL MODEL 

For the purposes of this study, the Sun, Mercury, Venus, 
Mars, Jupiter, Saturn, Uranus, Neptune, and Pluto are modeled 
as particles. The Earth is modeled as a triaxial rigid body 



and the Moon as an asymmetric rigid body. The Sun,. Moon, and 
planets interact gravitationally with translational and rotational 
motions fully coupled. Non-rigid, dissipative, and relativistic 
effects are not considered here but anticipated for future 
inclusion in the model. 


MATHEMATICAL MODEL 


A system of 41 second-order ordinary differential equations 
has been derived to represent the physical model described in the 
previous section. These may be summarized as follows: 


A. Motion of Sun with respect to center of mass of 
Solar System 


• • 


rj. = G 


IT 




( 1 ) 


B. Motion of Moon and planets with respect to the Sun 


r ± + Gtirq + m i ) 


r . 

1 


r . 


il 


11 

T 

j=2 


m. $ . u. . 

3 3 


(i = 2,3, ...11) (2) 


C. Rotational motion of the Earth 


{'e '} 




+ 


1 

2 


[3*3 



(3) 
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D. Rotational motion of the Moon 



The "short-hand” notation sip = sin($)' and ccf> = cos ( 41 ) 
has been utilized in the above equation. 

A full discussion of the above equations including rationale 
and derivations appears later in this report. 

Equations (1) represent the motion of the Sun with respect 
to the mass center of the Solar System as forced by the Moon 
and planets. 

Equations (2) represent the motion of the Moon and planets 

with respect to the Sun. The gradient is to be taken with respect 

to the coordinates locating each body with respect to the Sun. 

The force function U. . may be written as 

13 

U. . = U?. + U 1 (5) 

13 13 

P 

where LL ^ represents the mutual gravitational interaction between 
all masses treated as particles and U 1 is the mutual potential 
of the Earth and Moon regarded as finite rigid .bodies. 
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Equations (3) represent the rotational motion of the Earth 
with respect to a defined "reference" coordinate frame. The 
Euler parameters {B'} T ={3 q B 2 , 63 } represent the angular 

f 

deviations of the Earth from this "reference" frame. The forcing 
torques enter through the angular acceleration terms, 

0 

0)2 



The function (f(t)} defines the reference frame. A solution of 
these equations represents the deviation of the Earth from a uni- 
form rotation about a fixed axis in space, i.e., precession and 
nutation. 


Equations (4) represent the rotational motion of the Moon 
with respect to a defined "reference" coordinate frame. The 
Euler parameters {g ,,, } T = {go''/ 3 i * * ^ g 2.' ' • 3 3 * ' 3" represent 
the angular deviations of the earth from this "reference" frame. 
The forcing torques enter through the angular acceleration terms. 


0 

“l 

“2 

S3 


The rotation matrix [c(g ,f, )I and the terms in -AsiJ>, <f>, 

• ^ 

Actj) and their derivatives account for the motion of the "reference" 
frame. A solution of these equations represents the optical plus 
the physical librations. 


Rationale for Development of Equations 

The various existing analytical theories for the translational 
and rotational motion of the Earth and Moon are being looked at 
more and more critically . due to ever increasing observational 
accuracy (ref. 1) . 
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Many previously neglected effects must now be included in 
the mathematical models used to reduce the observational data. 
This, of course, leads to a better knowledge of these small 
effects and their causes. Currently unknown effects may also 
be discovered as the observational data is analyzed. 

The possibility of using LURE data to determine various 
geophysical and selenophysical parameters was pointed out in 
the introduction and in reference 2 . 

Accordingly, this report describes work undertaken to develop 
a mathematical model of the motion of the Earth-Moon system that 
has as few restrictions on accuracy as possible. Thus, the 
coupled rotational- translational motions of both the Earth and 
Moon are included in this model. 

A more immediate and specialized goal, however, is to be able 
to solve for the coupled rotational- translational motion of- the 
Moon for use in the reduction of LURE data to estimate the low- 
order gravitational harmonic coefficients of the Moon. Thus, 
this problem will be emphasized in this report. 

Several recent papers (refs, l, 3 , 4 , 5) have pointed out 
the facts that ( 1 ) analytical theories for the lunar translational 
motion and ( 2 ) analytical theories for the lunar rotational motion 
are not accurate enough to be used in the reduction of LURE-type 
data. Attempts are therefore being made to numerically integrate 
( 1 ) the equations of motion for the lunar orbit (ref. 3 ) , and ( 2 ) 
the equations governing lunar rotation (refs. 5,6). 

The above facts and the success of Cohen and Oesterwinter 
(ref. 7) in numerically integrating the motion of the solar 
system have prompted this attempt at a numerical integration of 
the equations of motion representing the coupled translational- 
rotational motions of the Earth and Moon. 

The formulation of both the translational and rotational 
equations of motion as a system of second-order differential 
equations was dictated by the general observation that Class II 
(second-order) numerical integration methods are more efficient 
(ref. 8 ) . 
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Since the force and torque evaluations at each integration 
step are very costly in computer time, it was decided to utilize 
Euler parameters rather than Euler angles in the rotational equa- 
tions. The relation of the rates of change of those parameters 
to the angular velocity components is algebraic rather than 
trigonometric in the case of Euler angle rates. Although two 
additional second-order equations are thereby added to the system, 
no trigonometric functions need be evaluated at each step. A time 
saving is thereby accomplished in the integration of the rotational 
equations. This approach has been common practice in the simula- 
tion of aircraft and gyroscopic motions (refs. 9,10). Advantages 
arise also in problem formulation and parameter estimation when 
Euler parameters are utilized (ref. 11). 

The reference axes used in the rotational equations were 
chosen so that the large angular rotation rates of the Earth and 
Moon with respect to inertial space did not have to be integrated. 
The reference axis for the Earth spins with respect to an inertial 
system about a fixed axis with a fixed rate equal to the mean 
sidereal rotation rate of the Earth. The reference axis for lunar 
rotation is centered at the Moon's mass center and its primary 
axis points to the Earth's center of mass. The axes of this system 
are parallel to' the unit vectors of a spherical polar coordinate 
system that locates the Moon with respect to a mean equator and 
equinox of 1950.0 rectangular system centered at tlie Earth. This 
approach is similar in philosophy to the Enke method of celestial 
mechanics . 

Coordinate Systems 

The coordinate systems utilized in this study are standard 
and are summarized in table 1 and illustrated in figures 1 and 2. 
Transformation between coordinate systems is accomplished using 
orthogonal rotation matrices in the sense 

{x'} = I R XX « 3 {x} 

where [R] is a 3 x 3 rotation matrix for a rotation of {x} 
into {x 1 }. 
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Table 1. Coordinate systems 




Axis 




Unit Vector 



Origin of 

Notation 

Fundamental 

Fundamental 

Secondary 

Notation 


No. • 

Frame 

(i - 1,2,3) 

Plane 

Direction 

Direction 

(i = 1,2,3) 

Remarks 

1 

Barycenter 

i g 

x i 

Ecliptic of 

Intersec- 

xy points 

j i # 

x 

Primary 


of Solar 


1950.0 

tion of 

toward 

inertial 


System 



ecliptic of 

North Pole 


reference 




1950.0 and 

of eclip- 


frame 





mean equa- 

tic of 







tor of 
Earth of 
1950.0 ' 

1950.0 



2 

Barycenter 

x. 

Mean equa- 

Same as 1 

X 3 points 

1. 

Secondary 


of Solar 

•L 

tor of 


toward 

- X 

inertial 


System 


Earth of 


North Pole 


reference 




1950.0 


of rota- 
tion of 


frame 






Earth of 
1950.0 



3 

Center of 

g 

x i 

Same as 2 

Same as 2 

Same as 2 

f! 

X 

Trans- 


mass of 




lating 


Sun 






frame ■ with 








respect to 
X. 

X 

4 

Center of 

Y i 

Same as 2 

Same as 2 

Same as 2 

5. 

X 

Reference 


mass of 




frame for 


Earth 






Earth rota- 
tion . 

Rotates at 
uniform 
rate of 
a with 
'respect 

H 

H 






(cont ’d. ) 

to x!. 




Table 1. Coordinate systems (concluded) 



Axis 



\ 

Unit Vector 


Origin of 

Notation 

Fundamental 

Fundamental 

Secondary 

Notation 


Frame 

(i = 1,2,3) 

Plane 

Direction 

Direction 

(i » 1,2,3) 

Remarks 

Center of 

, y i 

Plane of 

Axis of 

Axis of 

3 i 

Earth 

mass of 

equatorial 

minimum 

minimum 

"body 

Earth 


principal 

principal 

principal 


fixed" 



axes 

moment of 

moment of 


frame 




inertia 

inertia 



Center of 

Z. 

l 

Plane 

Axis points 

Axis points 

K. 

l 

Reference 

mass of 

formed by 

from Moon ' s 

opposite to 

frame for 

Moon 


radial and 

mass center 

longitud- 


lunar - 



longitud- 

to Earth ' s 

inal unit 


rotation. 



inal unit 

mass center 

vector as 


An "orbi- 



vectors of 


described 


tal ref- 



spherical 


in column 


erence 



polar coor- 
dinate sys- 
tem locating 
the Moon 
with res- 
pect to 
T i 

cen- 
tered at 
the Earth 
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frame" . 


Center of 

z . 

l 

Plane of 

Axis of 

Axis of 

it. 

Moon 

mass of 
Moon 

equatorial 

principal 

axes 

minimum 
principal 
moment of 
inertia 

minimum 
principal 
moment of 
inertia 

l 

"body 

fixed" 

frame 


The superscript T on a set of axes indicates a frame translating with respect to the unsuper- 
scripted frame. 


H 

to 







Equations of Motion 

Translational Equations , Reference 13 provides the equations 
of motion for the particles and centers of mass with respect to 
an intertial reference frame, viz. 


m. r. = 1 U 


(i =* 1,11) 


( 6 ) 


where 


U = G 


11 

ZE 

i> j=l 


m.m. 

*• . J. 
r ij 


(7) 


and 


mi 

is 

the Sun’s mass 

m 2 

is 

Mercury's mass 

m 3 

is 

Venus * mass 


is 

Earth's mass 


is 

the Moon 1 s mass 

m 6 

is 

Mars ' mass 

Itly 

is 

Jupiter's mass 

mg 

is 

Saturn's mass 

m 9 

is 

Uranus ' mass 

m l 0 

is 

Neptune's mass 

m X l 

is 

Pluto's mass. 


In equation (6) , 

?i = x ljL f, + x 2 . I 2 f x 3 . f 3 


(8) 


V i = ^ 3X X . + 3X 2i 


+ Ir 


3X 


3i 


(9) 


where X 


li, A 2i ' 


X^^ are the coordinates of mass 
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If the origin of coordinates is now translated to the Sun, 
the equations of motion for the Moon and planets are 








( 10 ) 


where 



The terms on the right-hand side of equation (10) arise from the 


force function U?j of equation (5) . Reference 12 provides 
U?. as 



since 



] ii] 
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The equations of motion for the Sun are 


• * XI 

- G Dm 
j=2 3 

where 

r. = X* . l[ + x' I2 + X* $3 
3 ID 1 2 J z 3] 6 

These equations follow directly from equations (6) and ( 7 ) . 

As will be discussed later, the mutual gravitational potential 
of the Earth and Moon, treated as finite rigid bodies, may be 
important to LURE accuracy. Thus, U in equations (6) and ( 7 ) 
should be of the form 

1111 m.m. 

u = G XX r 1 

i>j-l r ij 



Gm £} m 5 


y ± p ( s $) 

m=0 (r 45 ) 


• X cos ml ■+ Y sin mX 
nm nm 


Thus, if the origin is shifted to the sun, the resulting potential 
may be written as 



P I 

U. . + U 45 
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where 



j = 2 , r 11 

i / j 



/ G 
^45 




nra 


(S(j)) 


X 

nm 


cos mX + Y 


nm 


sin mX 


■] 




) ( 14 ) 


) 


Rotational Equations for the Earth . The rotational motion 
of a rigid earth must satisfy Euler's principal axis equations/ 
viz. 


( • 

“1 


/ \ 
M^A 


( 

kia) 2 w 3 

“2 

> ■ ( 

m 2 /b 

> - < 

k 2 a>iw 3 



m 3 /c 
\ > 


_k 3 ai 1 to 2 
V / 


where 

kj = (C - B)/A 

k 2 = (A - C)/B 

k 3 = (B - A )/C . 


In the above/ ok are inertial angular velocity components 
in the y^ frames. The moment components likewise are in 

this frame. K, B f C are the principal moments of inertia of 
the Earth. In order to orient the Earth with respect to the 
"reference" axes and the inertial axes the following sets of 
Euler parameters are introduced: 
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1) 


represents a notation from 


{ 3 } 

2 ) {£’} 
3) {3*'} 


represents a notation from 


represents a notation from 


{x!} 

to 


X 


{Y.} 

l 

to 


X 

'-r-' 

to 

(Vi) 


r 


, and 


F 


where 

(3) ~ L 3 0 , 3i7 32/ 3 3 J 


Reference 13 provides the following relation between the sets of 


Euler parameters 

representing the 

above successive rotations: 


{3 ’ 1 > = [3’] 

{ 3 } = 

[31 {3* 

} ‘ 


(16) 

where 

sS 


-32 

“33 




^ 1 

^ t 

T 

f 




Bi 

3o 

$3 

"02 



[3’1 = 

_ 1 

1 

r 

1 


(17) 


02 

-33 

0o 

• 3l ! 





1 

t 

I 




0 3 

02 

"3i 

3o 







J 



and where 








0o. 

"3i • 

-3 2 

-33 



[3] = 

3i 

3o 

“33 . 

32 


(18) 


32 

33 

3o 

-3i 




3 3 

”32 

3i 

3o 




The rotation matrices linking the above coordinate systems 

are 
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{ V - 

m x . Y ) 

{xp = 

[C O) 3 

(xp , 

(19) 

{y i } = 

'V 

- 

[c (8 * ) ] 

{Y i > , and 

( 20 ) 

(yp - 

IE X'y ] 

{x i - 

Ic(3") 

1 {x[} 

( 21 ) 


The rotation matrices introduced above have the following 
form where expressed in terms of Euler parameters {ref. 13) : 

2(6163 ~ 6 0 $ 2 ) 

2 ( 62^3 + 6061 ) 

2 2 0 2 
6 0 - 81 - 62 + 63 

( 22 ) 

The matrices [6' 3, ■[I 1 ], [c(g) 3 , etc. are all orthogonal 

and hence their inverses are their transposes (ref. 13) „ 

The angular velocity of the Earth can now be expressed in 
terms of the Euler parameters and^rates. The inertial angular 
velocity of the Earth is 

“ = “i3i + W 2 J 2 + w 3 3 3 (23) 


z z o 2 

6 0 + 81 - el - 83 
tc(6)] = i 2(6162 - 6063) 
2(8163 + 8082 ) 


2(8182 + 6063) 

2222 
60 - 61 + 82 - 63 

2 (6283 - 6061 ) 


or 


U W Y/X' + W y/Y 0 


(24) 


The reference axes, Y. , rotate at a uniform rate a about the 
J 3 axis so that 


W Y/X' = « J 3 


(25) 
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Reference 13 provides the following relation between the Euler 
parameters, their rates, and the angular velocity components 
in the rotating system: 


( ° ^ 

Wi 

-Y/X' 

( w 2 > = 2 tB * 3 (§’> (26) 

Y/X' 

0) 3 

Y/X' 

\ 

where 

{6'} T = LB 0 3i323 3 J 


To provide a consistent four-parameter representation of the 
angular velocity of the Earth, the vector “y/X 1 roust be projected 
on the {y^} axes and certain "augmented" matrices must be intro- 
duced. Accordingly, the angular velocity of the Earth assumes 
the form 


" 0 1 


' 0 1 

0)1 


0 


> = [c($')] a ( 

0 

> w 3 ; 




+ 2[g»] _1 


{$'} 


(27) 


where the "augmented" rotation matrix [c(3’)] a is of the form 



0 

_!jl 


0 

[c(8')] a = 

0 

0 

1 

1 

1 

1 

[c (8 ' ) ] 



0 

j 




(28) 


21 




where a = a + at. The combination 


[c(BM ] 


is now defined as {f(t)} and has the form 

0 

2 a 63 - 3o3i) 

2a (3 2$ 3 + Bq3i) 

a(B6 2 - el 2 - B 2 2 + B 3 2 ) 

A set of second-order equations can now be formed by differ- 
entiating equations (27) and solving for {g 1 }. These equations 
assume the form 


(f (t)}j 

(31) 


0 

o 

“1 

w 2 

U > 3 
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In equation 


( 31 ), 


7 0 ) 

0)i 

W 2 

0) 3 


{f (t)}. 


and 


r 0 ' 

0)i 

• > 
0)2 1 

0) 3 


are provided in equations (27) , (30) , and (15) respectively, 
term {f(.t)} is of the form 


The 


i * » 




{f(t)> = ( 


2a(BiB 3 + 81B3 ~ Bo$2 “ 0002) 
2 a(B 2 3 3 + B2B3 + B0B1 + B0B1) 
2 a (80 Bo — 81 Bi ~ B2B0 B3B3) 


and [B 1 ] is of the form 


[B * 3 = 


Vo 

- 0 i 

-02 

-03 

• I 

• t 

• T 

• ■ 

01 

Bo 

“S 3 

02 

• t 

• t 

* I 

• V 

£2 

03 

3 o 


b 3 

"02 

0 i 

• 1 

0 o 


) 


( 32 ) 


(33) 


Rotational Equations for the Moon . The derivation of the 
equations of rotational motion for the .Moon proceeds along a path 
similar to that used for the Earth, 

The rotational motion of a rigid Moon must satisfy Euler's 
equations, viz. 
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( 34 ) 


s • 

w « 




f ' 

M 1 /A' 


f . . . • 

kj W2 ‘W3 

• f 1 

0) 2 

> - 1 

M2/B ' 

> - i 

.1 1 » 

k2 Wi • W3 / 

• • 

0) 3 

J 


M3/C' 

i 


. » 1 I 

k 3 W2 j 


where all quantities here have analogous definitions to those of 
equation (15) . Primes are used to distinguish variables related 
to the Moon from those related to the Earth (un-primed) . 

The inertia rations k!^ in equations (32) have a more familiar 
notation, viz. 

k{ = a 

k£ » -3 (35) 

k 3 “ Y 


These ratios are related by the constraint 


ct — 


*Lz _ X. 

1 - 3y 


(36) 


Euler parameters are now introduced to orient the Moon with 
respect to its "reference" axes and the inertial frame. For this 
purpose, define 


{& * ' ’ } = { 


Dill 

3 ,ll 


I I I 


Dill 

3 3 


which represents a rotation from {Z.} to {z.}. The corresponding 

X X 

rotation matrix is 
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(37) 


{z i } = [R Zz ] {Z i } = ^ z i> - 


The inertial angular velocity of the Moon is 


to' 


U Z/X' + W z/Z 


(38) 


In terms of the Euler parameters { 3 * * 1 > and their' rates, the 
components of are 


f ° ' 

W lz/Z 

^ W 2z/Z ^ 

W 3z/Z 
\ i 


= 2 [B ' ' ' j 1 




(39) 


The angular velocity u z/x , of the "reference” frame is defined 
completely by the translational motion of the Moon with respect 
to the Earth* 


introduce the 


To determine the angular velocity to , , 

Z/ X 

spherical polar coordinates r, X, <p as illustrated in figure 
2. These are the coordinates of the Moon's center of mass with 
respect to 


{X.'J. 


Now, define the "relative" rectangular coordinates. 


V 


as 


Al 

= X 1 5 

- X 1 4 = 

r 

cos 

<P 

cos 

A 


a 2 

X 25 

- X 24 = 

r 

cos 

4 > 

sin 

X 

(40) 

A 3 

X 3 5 

~ X 3 4 “ 

r 

sin 

<t> 

0 
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Inverting these expressions provides 


r = + A 2 + A3 


= tan^ 1 (A 2 /A x ) , 


tan" 1 (L^/Ai + A 2 ) 


Since the unit vectors k^ are related to those of the spherical 
polar system by 


the inertial angular velocity of the axes {Z^} can fc> e written as 


“z/X* “ ^ 4> = ^ ~ ^ 


The above vector may be projected on the {Z^} axes providing 


. - X [, 


-V -y *1 # 

cos (j> k 3 - sin <j> k x + <J> k2 


The components 


= -A sin <|» 


= A cos <J> 
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can be related to the relative position coordinates A . and the 
relative velocity components A^. To do this, differentiate 
equations (40) , solve for A^, and put the results in the 
matrix form 


• 


c$cA 

-sA -cAsiJ) 


( • 
r 

• 2 

M 

> = 

ccj>sA 

S(f> 

cA -sAscj) 

0 ctj) 

( 

rAc^ ) 
; r * , 


(45) 


Equation (45) can be inverted to provide 



f . 


p 


- 


\ 

• 


r 


c<f)cX 

sAc<}> 

C(j) 


Ai 


• 

rActj) 

* - 

-sA 

cA 

0 

\ 

a 2 > 


r 4> 


-cAs<f) 

-sAs<j> 

ctf> 


£3 


\ i 






\ / 

In the above equations, the 

"short- 

■hand" 

notation 

and s<J> 

- sin 

c p has 

been utilized. 




(46) 


Now, the components of ca„ /v ., in the { 2 .} frame are 

lx/ X 1 


w iz/x' 




r • 

( w 2Z/X ' 

> 

= [c(3'")3 


- 1 > 

[ U 3Z/X ' J 




* 

Ac<}) 


(47) 


• • 

The quantities -As^, <£, Ac< j> follow from equations (46) and 

(40) as follows: 
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£~ A ! cAs<f) - A 2 sAs<}> + A 3 c<{>] 


( 48 ) 



[-AjSA + A 2 cA1 


where 


sAc<|> = A 2 /r 
s4> = A 3 /r 
cAc<|) — Aj/r 


c<j> 

sA 



The absolute angular velocity of the Moon can now be obtained 
by adding (39) and (47) to obtain in augmented form. 



Equation (49) can be solved for {p’*’} and differentiated 
to obtain the second order equation for {p 1 ''}. Thus, 
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In equation (50), [ 3 ’ * 1 ] is of the form given in equation 

(18); [c( 3 ,,, ) 3 a the f° rn * given in equations ( 22 ) and (28); 

[g* ' ' ] is of the form given in equation (33); 10 w{ WgJ are 

given in equations (34); L0 co{ w 2 w|J are given in equations (49); 
and the elements of [c(g ,,, ) 3 A are 


Cll A Cl2 A Cl3 A Ci4 A ° 21 A Csi A ° hl A 0 


c 22 a ~ 2 (8 0 8 0 + B]_8i - 8282 ~ 33^3) 

° 33 a = 2( Mo “ 3 1 3 1 + 8282 “ 8383) 

= 2(3 0 3 0 “ 8181 - 8282 + 8383) 

C 4 4 A 

c 23 a ~ 2(8 1 3 2 + 8182 + 8083 + 8083 ) 


(51) 

(cont 'd. ) 
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(51) 
(eoncl *d. ) 


c 24 a = 2 ( 3 1 3 3 + Pi33 “ ” ^0^2^ 

c 3 2 A = 2 (3 i32 + 3 1 3 2 ~ 3 0 3 3 “ 3 0 3 3 ) 

c 3 4 a = 2(g 2 33 + 32?3 + 3o3l + 3 0 3 1 ) 

C 4 2 a = 2(3^3 + 3i33 + 3 0 3 2 + 3 0 3 2 ) 

Clf3 A ~ 2 (3 2 3s + 3203 - 3o3i ** 3o3i) 

All g's- in the above equations have a triple-prime superscript. 

The derivatives of the polar coordinates and their rates in 
equation (50) are 


& (Ac*) =Ai.[-Sii + ki] + A 2 [-S|i- 


V sX , cA 
Al + a 2 


ft “ 

rcX 
r^ - 


(52) 


dt 


(d) = A I + rcAstj) $c\C(j) I 

9 1 L r r 2 ~ r J 


+ a 2 - * 


AcXstj) + (j>s\c(p + rs 
r r ' 


rsAstj) j 

r2 J 




r • # n 

(53) 

+ A, _ _ rc£ 


3 L r J 


v cXs* _ v sXs<t + v c£ 

1 3T ^ ^ 


M - if {£ <***>} • 

(54) 
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where 


A. 

i 




• • 

A. 

x 




( 55 ) 


are available from the integration of equations (10) . 

Forces and Torques 

The gravitational forces between a set of particles was given 
in equations (10) . If the Earth and Moon are considered as rigid 
bodies then mutual gravitational forces and torques arise and 
must be modeled. Also torques exerted on the Earth and Moon 
due to a point mass Sun must also be considered. These are 
developed in this section. 

Mutual Gravitational Potential . There are two approaches in 
the literature for deriving the mutual forces and torques. Approach 
(A) involves deriving a mutual gravitational potential and then 
finding the gradient of this potential with respect to the trans- 
lational and rotational variables to give the forces and torques 
(refs. 14, 15). Approach (B) involves a direct integration of 
a differential force and torque over both bodies (ref. 16) . 

Approach A appears to be more easily developed when higher 

order gravity harmonics' than the second are included for each 

body. Approach B is more concise than A for the case when only 

second order terms are retained for either one or the other body. 

* 

Also, the effect of Earth oblateness on lunar torques can readily 
be derived using this approach. 

For generality and ease of extension to higher orders, 

Approach A will be- followed here. The concise results of Approach 
B are presented in Appendix A. 

For the purposes of this report the Earth and Moon will be 
modeled as follows: 
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Earth (ref. 17, {y. }) 

l 

c 20 = -1.082637 x 10- 3 

c 2 1 = s 2 1 = 0 
c 2 2 = 1. 5362 x 10~ 5 
s 22 = -8.8149 x 10~ 7 

These' values provide the following moments of inertia: 

A = .33912 Ma 2 
B = .33906 Ma 2 
C = .34017 Ma 2 

if the dynamical flattening H = (C - A) /c = 3.27293 x 10 3 
is adopted from reference 20. 




Moon (refs. 

18, 19, 

{z ± }) 

c 20 

= 

-2.0272 x 10 - 1 * 

(ref. 

18) 

C 2 1 

= 

s 2 1 = 0 

(ref. 

19) 

c 22 

= 

2.221 x 10“ 5 

(ref. 

18) 

s 22 

= 

0.0 

(ref. 

19) 

c 3 0 

= 

3.9 x 10- 6 

(ref. 

19) 

^31 

= 

28.6 x 10“ 6 

(ref. 

19) 

c 32 

= 

6.0 x 10~ 6 

(ref. 

19) 

c 33 

= 

2.7 x 10~ 6 

(ref. 

19) 

s 3 1 

=r 

8.8 x 10“ 6 

(ref. 

19) 

S32 

= 

1.8 x 10~ 6 

(ref. 

19) 

s 3 3 

=5 

-1.4 x 10- 6 

(ref. 

19) 

O 

& 

0 


23.3 x 10~ 6 

(ref. 

19) 
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C 4 1 

= 

11.1 X ' 

10“ 6 

(ref. 

19) 

c 42 

= 

-2.48 

X 

10“ 6 

(ref. 

19) 

c 43 

- 

-0.17 

X 

10“ 6 

(ref. 

19) 

C 44 


-0.25 

X 

10~ 6 

(ref. 

19) 

s 4i 

= 

-2.61 

X 

10~ 6 

(ref. 

19) 

s 42 

= 

-3.28 

X 

10“ 6 

(ref. 

19) 

s 4 3 

= 

-0.45 

X 

10~ 6 

(ref. 

19) 

s 44 

= 

0.27 : 

X 

10“ 6 

(ref. 

19) 


These values provide the following moments of inertia: 

A f = .391753 M'a' 2 
B* = .391842 M*a ' 2 
C' = 0.392 M'a ' 2 (ref. 18) 

if the values of 8 and y are taken to be . 

8 = 631.1 x 10" 6 
y = 226.8 x 10-6 


as in reference 18. 

Reference 15 provides the general form of the mutual potential 
between two arbitrarily shaped rigid bodies in the form 


00 ' 

u Is ~ E Up = G'r " 1 
£= 2 * 


CO 


E 

n=2 


n 

E„ 

m=0 



P 


nm 


(sin <J>) 


[X cos ml + Y sin mX] 

L vim -v^m 


nm 


nm 


}■ 


(56) 


33 



0 


... . p I 

where the term Gr 1 has been included xn and where Uj = 

due to the choice of coordinate system. Here, M is the mass of the 

Earth and M* is the mass of the Moon. The X and Y are 

P a nm nm 

functions of a. a*, c , c' , s , s' where a and a' 

’ ' pr qs pr' qs 

are the mean equatorial radii for masses M and M* and the c's 

and s 's represent the harmonic coefficients for both bodies. 

Reference 15 provides the lower order values of X and Y as 

nm nm 

follows: 


X 2j = a2c 2j + a * 2 °2 j 


(j = 0,1,2) 


Y 2j = a S 2j + a ’ S 2j (j = X ' 2) 
X 4 q = 6a 2 a* 2 (c 20 c 21 “ ^C2x t - : 2i — 


+ 2c 2 2 c 22 + 2 s 2 2 s 22^ 


X41. = 3 a 2 a ,2 (c 21 c^ 1 + c 21 c£ 0 - c 21 c£ 


“ C 22 C 2 x + S21®22 ~ a 22®21^ 


Y 41 = 3 a 2 a* 2 (c 20 S22. + s 21 c 20 " c 21 s 2 


- S 22 C 21 + S 21 C 22 + c 22 s 2l) 


X 42 = a 2 a" 2 (c 20 c 22 + c 22 c£ 0 + c^c^ - s 21 s^) 


Y 42 = a 2 a ' 2 (c 2 qS 22 + s 22 c£ 0 + c 21 s£x + s 21 c£i) 


“ \ a 2 a' 2 (c 21 c 22 + c 22 c^ - s 21 s£ 2 - s 22 s^ L ) 
= — a 2 a ,2 (c 21 s 22 + c 22 s 21 + s 21 c£ 2 + s 22 c 2l) 


— a 2 a' 2 (c 22 c 22 - - s 22 s 22 ) 


= k a 2 a 1 2 (c 22 s 22 + s 22 c 22 ) " 
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The coordinate systems implicit in equation - (56) is illustrated in 
figure 3. 

If {y.) is chosen as the primary reference, then c , s 

are the harmonic coefficients of the Earth referenced to {y . } 

T * 

and c* , s' are those of the Moon referenced to {y. }. Like- 
qs qs . l 

wise, r, X, and <p are the spherical polar coordinates of the 
lunar mass center with respect to {y.}. If the {z . } axes are 
chosen as the primary reference, then the above quantities are 

m 

referred to {z^} and {z^ }„ 

The relative orientation of the above axis systems can be 
expressed as follows: 


{yp 


a a' a’ ' 

3 3' 3” 

Y y ^ y ^ ^ 


{ z^ } = 


(58) 


with a more precise definition of the a's, 8's, and y's to 
be given later. 

Force on Moon Due to Earth. In terms of spherical polar coor- 
dinates r, A, <f> locating the Moon with respect to {z^ }, the 
vector force may be written as 


£1 


m 5 


-j. ^ 9^45 ^ ^ U 4-5 , f A 

.3r~ ^r + r cos 4> dX X \ + r 3<j> ^ 


(59) 


The general mutual potential may now be specialized to the problem 
at hand as follows: 
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( 60 ) 


■Now, 


U 9 = Gr 


U 3 


|^(a/r) 2 jc 20 P 20 ( s< f>) + p 2 i(^ 2 i c ^ + s 21 sX) 

+ P 22 (C22 C2?1 + s 22 s 2X) I + (a*/r) 2 jc£ 0 P 20 

+ P 22 ( c 22 c2 ^ + S 2 2 s2 ^)| J 

= Gr - 1 j\a'/r ) 3 |c^ 0 P 30 + P 31 (chcX + s| lS X). 

+ P 32 {c 3 2c2 X + s| 2 s2X) + P 33 (c$ 3 c3X + s^ 3 s3X)|J 


U k = Gr 


-1 


(61) 


£ (a'/r) 4 ■|c 4 q p 40 + P41 (C41CX + s^isX) 

+ P42 (c42 c 2X + S4 2 s2X) + P43 (c4 3 c3X + s 4 3 s 3 X) 

+ P44 (C44C4X + S44S4X) j- + (a 2 a ,2 /r 4 ) j 6P4 0 (2c 22 C22 
+ 2 s 22 S 22 ) + 3P 4 1 ( ^2 1^2 0 ~ c 21 c 22 + S 21 s 22 ) c * 

+ 3 P 41 (s 2 iC^ 0 - C 2 1 S 2 2 + S2iC^ 2 )sX 

+ P42( c 20 c 22 + c 22 c 2 0 ) c 2 -X + P42( c 20 s 22 + s 22 c 20 ^ s2 ^- 

+ (P 43 /2) (c 2i c^ 2 “ s 2 iSi 2 )c3X + (P 43 /2(c 21 s£ 2 

+ s 2 1 ^22 ) < -' 3 X + (P44/2) {C22C22 — S 22®22^ C 4 X 
+ (P44/2). (C 22 S^ 2 + S22C22) S 4 X |J » 


(62) 


.the partial derivatives of U 45 can be calculated and are 


3U 2 /3r = -3G a 2 r ~ 4 { c 20 P 20 + P 2 l (c 2 icX + s 21 sX) 

+ P 2 2 (c 22 c2 * + s 2 2S2X) I - 3G a ’ 2 r“ 4 jc£ 0 P 20 (63) 

+ P 2 2 (c 22 c 2X + S 22 s2X)| 

3U 2 /3X = G a 2 r ~ 3 jc 20 P 20 + p 2 l (“ c 2 l s ^ + s 21 cX) 

+ 2 P 2 2 (“C 2 2 sA + s 22 c2X) [ + G a , 2 r ~ 3 |c£ 0 P 20 (64) 

+ 2P 22 (-c 22 s 2X + s£ 2 c2X)V 
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3u 2 /3<j> 


= G 
+ p 
+ p 


a 2 r ~ 3 {c 20 P 20<|) + P 21(J> 
22 ^ *6 S 22 s2A) ^ 

22<p ( c 22. c 2^ + s£ 2 s2X)j- 


(c 21 cX + s 21 sX) 

+ G a 2 r ~ 3 |c£ 0 P 20< £ 


3U^/3r — -4 G a’ 3 r~5 |c^ 0 P 30 + P 31 (c| 1 cX + s| lS X) 

+ P 32 (c 32 c2X + s 32 s2X) + P 33 (c 33 c 3X + s 33 s3X)| 

3U 3 /3X = G a ,3 r~ 4 |c 3 qP 30 + p 3 i(~C 3 ].sX + c 31 cX) 

+ 2P 32 (~c 32 s2X + s 32 c2X) 

+ 3P 3 3 (-c 3 3 s 3X + s 3 3 c3X) | 


3U 3 /3<j) = G a 1 3 r ~ 4 {c$ 0 P 30(j) + P 3 x^ ) (c| 1 cX + s^sX) 

+ ^*3 2 ^ ( c 3 2 c2 ^- + s| 2 s2X) + P 3 3 ^ (c 3 3 c3X + s 3 3 s3X) 


3U 4 /3r = -5G a , 4 r " 6 jc^oPtto +-P 41 (c 4 icX + s 41 sX) 

+ P 42 (c 4 2 c 2 X + s 4 2 s^X ) + P 43 (c 43 c3X + s 43 s3X) 

+ P 44 (c4 4 c4X + s 44 s4X) | - 5GMM'a 2 a , 2 r ~ 6 j 6 P 40 X 40 

+ 3P 41 (X 41 cX + ? 41 sX) + P 42 (X 42 c2X + Y 42 s2X) 

+ (P 43 /2) (X 43 c3X + Y 43 s3X) 

+ (P 44 /2) (X 44 c4X + Y 44 s4X)| 


3U 4 /3‘X 


= G a ,4 r~ 5 ^ c 4 qP 4 q 4- P 41 (-c 41 sX + s 41 cX) 

+ 2P 42 (-c 42 s2X + s 42 c2X) + 3P 43 (-c 43 s 3X + s 43 c3X) 

+ G a 2 a' 2 r“ 5 |6P 40 X 40 
+ 3P 41 (--X 41 sX + Y 41 cX) + 2P 42 (-X 42 s2X + Y 42 c2X) 

+ ,(3P 43 /2) (~X 43 s3X + Y 43 c3X) + 2P 44 (-X 44 s4X 
+ Y 44 c 4X) j. 


+ 4P 44 (-c 44 s 4X + s 44 c4X) | 


( 65 ) 


( 66 ) 


(67) 


( 68 ) 


(69) 


(70) 
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3uJ/9cf> = G r“ 5 | P lf0(J) (a ,4 c4 0 + 6a 2 a ,2 X 40 ) 
+ P 41 [(a ,4 c4i + 3a 2 a' 2 X 41 )cX 

+ (a ,2 si x + 3a 2 a ,2 Y 41 )sxJ 
+ P 4 2 ^ £{a ,4 c 42 + a 2 a ,2 X 42 )c2X 
+ (a' 4 s 42 + a 2 a* 

+ P 43 <J> 1+043 

+ {a ,4 s 43 + 2 Y 43 )s3xJ 

+ P 4 4 (j) [(a ,4 ci 4 + | X 44 )c4X 

+ (a ,4 s4 4 + j Y 44 )s4xj | „ 

In the above equations, 

X 4 q — 6a 2 a 1 2 [ 2022^22 2s 22 S 22 ] 

X 4 1 — 3 a 2 a* 2 [c 2X c 2 g “ c 2X c 22 "** ®2l^22^ 
Y 4i = 3 a 2 a' 2 [s 21 c£ 0 - c 2 is£ 2 + s 2i c 22l 
X 42 = a 2 a ' 2 [c 20 c 22 + c 22 c 2 o 3 
Y 4 2 “ a 2 a 1 2 [c 2 q s 2 2 + s 22 c 2 0 i 
x 43 = ~2 a 2 a ,2 [c 2 iC22 ~ S 21 S 2 2 1 
¥43 = J a 2 a ’ 2 [c 21 S22 + s 21 c 2 2 3 
x 44 — 2" a 2 a 1 2 [c 2 2 c 22 "* ^22^22^ 

Y 44 = ~ a 2 a 1 2 [c 22 S 22 + s 22 c| 2 ] 


2 Y 42 ) s2Xj 
+ i X 43 )c3X 
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and 


= 3si P c <t> 

P 21(j) « 3c (24.) 

P 22 (j, - -Sstpccp 
P 3 0 ^ = | c<j>(5s 2 <f, - l) 

P 31(j) = | S$ + S4>(2 c24, - S 2 (J>) 

P 32 ^ ~ 15c<j> (c 2 (j> — 2s 2 <{>) 

P 3 3 ^ = -45c 2 4>S(J) 

P 40(j) » ~ (140s 3 4) - 60s4>)c4> (73) 

P 4 ]_^ — [(c 2 , — s 2 <f>) (7 s 2 4> — 3) + 14s 2 4>c 2 4>^j 

P 42 ^ = ~ |jL 4 s 4> c 3 4> - 2s(f>c<}> ( 7s 2 <j) - 1)J 

P 43 = 105[c 4 <j> ~ 3s 2 <j)C 2 ^3 

P 44 ^ - -420c 3 <j>s<}> 


Since the harmonic coefficients c. . and s. . are referred 

T ■ 1J 1J 

to a coordinate system } that rotates with respect to the 

earth, they are functions of time. These functions may be evaluated 

by noting the definitions of the c. . and s. . (ref. 21) in terms 

r j r 3 

of the inertia integrals, viz. 
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c 2 0 ~ 

C21 = 
S 2 1 = 
c 22 ~ 
S 2 2 = 


1 

a 2 M 


1 

a 2 M 

1 

a 2 M 

1 

4a 2 M 
1 

2a 2 M 




I 11 + I 2 2 
_ 




( 74 ) 


and by noting the transformation laws of the inertia matrix, viz. 


{z.} 

II 

1 1 

H3 

{yp 





"1 - 
z 1 2 1 

I 

2 1 2 2 

Iz l z 3 


'a 

0 

0 ' 

I 

2 2 2 1 

. I 

Z 2 Z 2 

I 

z 2 Z 3 

II 

*=> 

UJ 

i-3 

0 

B 

0 

I 

L 2 3 2 1 

I 

z 3 Z 2 

I 

z 3 Z 3 _ 


0 

O 

c 


The above equations provide the desired functions as follows: 
a 2 Mc 2 o = -|[A + B + C-3(ct ,,2 A+ 3 " 2 B + y f,2 C)J 
a 2 Mc 2 i = aa* 'A + Bf^'B + YY ,f C 

a 2 Ms 21 = a'a' 'A + B^'B + y'Y ,,c (77) 

4a 2 Mc 22 — A(cj’ 2 - a 2 ) + B(3' 2 - 3 2 ) + C(y' 2 - Y 2 ) 


2a 2 Ms 22 = aa'A + 33'B + yy'C 
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'Finally the components of the force on the Moon due to the 
Earth must be found along the {X^} axes. Accordingly, 


5 x! 

x 


m. 


3U 


45 


+ r -5? 


3r 


3U Is 


(I 




1 9U ^ 
V + ^ -2 T- 






x i) 


( 78 ) 


Torque on Moon Due to Earth . This torque can be derived by 
first expressing the mutual gravitational potential U ^ 5 in terms 
of the direction cosines relating the rotational orientation of 
the moon to the earth; and then by calculating the moment components 
as follows (ref. 22 ) : 


tt-t M = ct 1 ’ U 1 , - a'U 1 , , + 3' 1 Uq, - 3' uj, , 
MM' z-i a a ' 3 3 


+ Y 1 U*. - y* 


rf^r M = a U 1 ,, - a'* U 1 + 3 U* , 
MM' zz a a 3 


- 3 


1 ' ul 


+ Y - Y’* U* 


(79) 


M = a 1 U 1 - a U 1 ' + 3 1 U? ‘ - 3 U? , + y ' U 1 - y U 1 , 

MM 1 z 3 a a 3 3 Y Y 


A derivation of the above relations is presented in Appendix C. 

X 

In order to derive the torques, the term U 2 and the second 
order coupling terms in U 4 [see eqs. (81) and (82)] will be 
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treated together. Finally the term U 3 and the remaining terms in 
U ©« will be treated [see eqs. (86) and (87)]. 

The reference axes for the second order and coupling terms 
are {y . } . Thus the c ! . and s ! . are functions of the orien- 
tation angles. These functions are 


-I* 


c 20 “ 


2 M* 


yiyi + I y2y2 _ 
2 


C y 3 y 3 _ 


,2 m' 


A + B + C 3 , , 0 1 . r-> 1 i2 , n i 

2 “ ~2 (A'y^ + B'y ,z + C'y 




cii 


1 

a ' 2 M' 


X yly3 


a' 2 M’ 


[ayA 1 + a'y'B* + a ,, y ,, C , l 


• 1 x * 
S 21 _ x 


. 2 M . y3y2 


(80) 


2 M' 


[yBA* + y * B *B ' + y ,, 3 ,, C'] 


C22 = 


4a 


— — [a. 
1 2 m* L 


(8 2 - a 2 ) + B* (8 1 2 - a' 2 ) 


C‘ (3* 12 - a* ,2 )j 


s 22 ~ 


4a 2 M 1 


[«8A’ + a'B'B' + «' '8* 'C’] 
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. I I 

The potential U 2 + COU pii n g thus assumes the form: 


U2 — G r ~ 1 j^(a/r) 2 {c 20 P 20 + c 22 P 22 c 2 X} 

+ (a'/r ) 2 { c 2 qP 2 0 + p 2i ( c 2lcX + s 21 sl) 


+ P 22 (c 22 c 2 X + s 22 s 2 A) 


coupling 


= G a 2 a ,2 r 5 ^P 40 | 6 c 20 c 21 + 12 c 22 c 22 ^ 

+ P 41 (|“ 3c 22 c £l} cX + |3c 20 sii + 3 c 22 s^i J- sX^ 
+ P 42 ^{ c 20 c 22 + c 2 2 C 2 0 } °2X + j c 20 s 22| s2X^ 

+ p 43({§- C 22 C 21 } c3 X + { ” c 22 s£x £ s3X^ (8: 

+ P 4 4 c 2 2 c 22 }c4X + C 22 s£ 2 ^ S4x)l 


r* cl* l c 4 "X 


Now, using equations (79) to (81) the torque components due 


to the terms in U 2 are; 


M 

zi 

= GMr~ 3 (B 1 - C') ["“3P 2 qY 1 Y ' 1 

+ P 21 cX (a'Y*' + a ,, Y l ) 

p2?^2X 


+ P 21 sX(8 , y” + Y'S") + 2 

(B 1 8* ' - a 'a 1 ' ) 


P 22 s2X "I 

+ (a * 8 1 1 + 8'a’ ') 

(83) 


r* - 

(cont'd. ) 

M 

z 2 

= GMr- 3 (C* - A 1 ) 3P 2 0 YY ' ' + 

P 22 c2X 

P 21 cX (aY 11 + YCi* 1 ) 


+ P 21 sX(8y ,t + y3") + 2 — 

P 2 2 s 2X "l 

+ 2 (aB 11 + Ba'*) 

(BB ' 1 ~ aa") 
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= GMr “ 3 (B 1 - A*) |^3P 20 YY + p 21 c ^(Y a ' + CXY*) 

P 2 2C2X 

- P 21 sA{y 3' + y'3) + 2 («“’ “ 

P 22 s2A -j 

- ^ (Ba* + aB’) 


(83) 
(concl- 1 d. ) 


The above components are for an arbitrary orientation of {y^} 
with respect to {z^}. In the derivations presented in the litera- 
ture, the Earth is treated as a particle so that the relative orien- 
tation of- {y i } and {z^} is immaterial. To recover those results 
a special orientation of the {y^} . may be taken. If y^ is taken 

yi 

to be pointing at the Moon then £ = c<j>cA = — =1, m = c<psX = ■ 

y 2 y 3 

~ = 0, and n = s<{> = — - = 0. Also, £* = a, m 1 = a * , and 

n f = a 1 ’, where Jl'/Ei^n 1 are the direction cosines of the Earth 
with respect to {z^}. This reduces equation (83) to a more 
recognizable form, viz. 

M x = 3GMr~ 3 (C' - B'Jm'n' 

M z2 = 3 GMr ~ 3 (A* - C')Jl'n' (84) 

M z3 = 3GM (B 1 - A'H'm' . 


Actually, equations (83) and (84) are identical as algebraic 
manipulation will show. 


The coupling terms in 



are handled similarly. 


Thus, 
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M 


Zl 


= GMr“ 5 a 2 (B* — C ' ) |^P 40 | c 20 (a* *y ' + y^a') 

+ c 2 2(3 , 3 ,? - a’a ,, )}+ 3Piti {-c 2 2(a ,, Y , + y^a'JcA 

+ sA (y ' 3 * 1 + y”3 l ) (c 20 + c 22 ) \ 

P 42 

+ c 20 | c2A (a'a 1 1 - 8 ' 3 1 f ) + s2A(a"3' + a'3")) 


M 


z 2 


P 43 

- 3P 42 c 22 c2Ay 1 y 1 ' + — c 22 {c3A(a ,, Y* + a'Y' 1 ) 

. P 44 

+ s3A (3 1 ' Y ' + 3 ' Y' * ) } 2 — c 22 |o4A(e , e" - a'a* 

+ s4A (a 1 '3' + 3' ’a')} J 

- GMr“ 5 a 2 (C 1 - A’) j^6P 40 {c 20 ( Y ,, a + ya") 

+ C22{33' ' - aa' ' ) } + 3P 41 | -c 22 (ay* 1 + ya'.'JcA 

P 42 

+ sA(c 20 + c 22 ) (By* ' + y3' ') } + ~ 2 ~ c 2 0 {(33" 

- aa ,, )c2A + (a3" + 3« ,, )s2A|- -3P 42 c 22 c2Ayy 1 * 

p 4 3 

+ — 2 ~ c 22 { ( a Y 1 f + ya'McSA + (3 y" + y3")s3A{ • 

+ c 22 {(33" - aa'')c4A + (a3" + 3a ,, )s4A}J 


M Z3 - GMr 5 a 2 (B 1 - A’) |^ 6 P 40 {-c 2 Q(a'y + y'a) 

+ c 22 (aa' - 33 ' ){ + 3 P 41 { + c 2 2 cA(a'y + y'a). 


42 


“ (c 2 o + C2 2 )sA(3'y + y'3)}. + — 2 ~ c 2 g |c2A(aa 
- 33') - s2A(a3* + 3a') } - 3P 42 c 22 c2Ayy' 

P 4 3 

+ -tt- c 22 {-{a'y + y'a)c3A - (3'y + y'3)s3A{ 
c 22 { + (aa 1 - 33')c4A - s4A(a'3 + a3')}J 


2 U22 
P 4 4 


( 85 ) 
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Finally, the torque components due to the U 3 potential 
X 

terms and the U 4 potential terms not already treated, viz. 


0 ? = 2 
3 j- 


{ p 3oW c 30 + p 3 iW(c^i cos X + S 3 X sin X) 

+ P 32 (<J>) (C 32 cos 2 A + s 3 2 sin 2 A) ( 86 ) 

+ P 33 (<J>) (c | 3 cos 3X + s § 3 sin 3X) 

U 4 = p | p 40 c 40 + P 4 1 (<J>) (C 41 cos X + S 4 1 sin X) : 

+ p 42 (<j>) (c 4 2 COS + s 4 2 s i- n 2X) 


+ p 4 3 ( 0 ) (c 43 cos 3X + si 3 sin 3X) 

+ P 44 (<{>) (c^ 4 cos 4A + s4 4 sin 4X ) } 


(87) 


will be derived. 

The general approach for calculating torques used earlier will 
now be modified. The Earth may be considered to be a particle in 
calculating the torques on the Moon arising from third and fourth 
degree terms in the lunar potential. Reference 23 provides a simple 
approach based on this fact that will be followed here. 

Consider that the reference axes in equations ( 86 ) and (87) 
are the {z^} axes. Thus the c |j ' s an ^ s ij* s are constant 
and r, <J>, X are the spherical polar coordinates of the center 
of mass of the Earth with respect to the axes {z^}. 

For a point mass Earth, the torque on the finite Moon due to 
the Earth is equal and opposite to the torque on the Earth due to 
the Moon. Thus, reference 23 finds the torques to be 
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M 

Z 1 


GM 

A 1 



3(^3 + uj) 

3z 2 



GM 
B 1 


■ Zl a(pj + 0») 

3z 3 


M 

Z 3 


GM 

C* 



3(Ug + O 4 ) 
3 z j 


z 2 


z 3 


z i 


3 (U 3 + U 4 ) ] 

3z 3 J 

3(U 3 + U l) 1 

3i7 J 

stu 1 + u 1 ) 1 

3z 2 J 


( 88 ) 


where 


\ 


f 


/ \ 

Zl 


cos <J> cos A 


£1 

z 2 

> ' * < 

cos <f> cos A 

> -= < 

*2 

z 3 

> 


sin <p 

i J 


*3 

\ / 


(89) 


Reference 23 then provides 



3 GMM'a' 3 r~ 4 T n ,, 

2 A 1 *2 (- 1 - 


5^1)c 3 o “ 10£ 1 £ 2 £ 3 c 3 j 


- 10* 2 (X •+* *| - 2&f ) c 3 2 - 60£ 1 £ 2 £ 3 c 33 


- £3 (1 - 5*| + 10A|)s 31 + 20£ 3 (£ 3 - £|) s 32 


(90) 
{ cont ' d . ) 


+ 30£ 3 (£ j - 


^ 2 .) s 33 


] 
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M_ 


3 GMM'a ,3 r- 4 
2 B' 


h 


(1 - 5 £ 3 )c 


3 * c 3 0 


+ £ 3 (X ~ 5£ 3 + X0£j)c 3 ^ 


- XOS-! (X + Jtf - 2£?)c 32 


- 30£ 3 (£ 1 - £2)^33 + X 0 £ 1 £ 2 £ 3 S 31 


+ 20 £ 2 (£j ~ £3)332 “ 60 £ 1 £ 2 £ 3 s 3 3 


+ - 


3GMM'a f 4 r“ 5 


B ! 


[-} 


£fc 41 4- 3 s£ 1 C4 3 


] 


(90) 
(concX’d. ) 


w 3 GMM'a ,3 r- 4 

M z = 2 — c* 


[- 


£ 2 (X " 5£ 3 )c 3 i 4- 40£j £2^3 c 32 


4- 30£ 2 (3£i " ^ 2^ c 33 & l ( 1 — 5£ 3 )s 3 j 


- 20£ 3 (£i - £|) s 32 - 30£ 1 {£j - 3£f)s 33 


C* 


3GMM'a ,4 r~ 5 f 2 4 1 

+ J — - — 5£ 1 s If2 - X40 £iS 44 


In summary, the force components in the inertiaX frame {X^} 
may be found from equations (78) and the preceding definitions found 
in equations (59) to (78). The torque components on the Moon- due 
to the Earth resoived aiong the {z^} axes are the sum of the 
torques in equations (84), (85), and (90). 

Torque on Moon Due to Sun . Since the Sun is treated as a 
particXe, the torque exerted on the Moon is of the same form as 
equations (84), viz.. 
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( 91 ) 


M = SGMjrup 3 (C* - B f ) 

m Q 

n © 

M = 3GM 1 r 14 “ 3 (A 1 - C') 

*6 

n o 

M = 3GM 1 r il T 3 (B' - A*) 

Z 3 1 




where £’ , ml . nb , are the direction cosines of the Sun with 
© 0 0 

respect to {z^}. 

Force on Earth Due to Moon . This force is derived in a manner 

analogous, to the force on the Moon due to the Earth presented 

T 

earlier. The reference axes are chosen to be the {y^ } set so 

that the c.. and s.. are constants and the c!. and sf. vary. 
ID rD X] X] 

An assumption made xn the following is that only second degree terms 
in the lunar potential are important to the motion of the Earth. 

• The mutual potential then assumes the form 


U 45 = Gr 


- 1 £ (a/r 1 ) 2 {c 20 P 20 (s4>) + ?zz (s<j>) (c 22 c2A 

+ s 22 s 2A)} + (a ’/r) 2 { c£ 0 P 20 (s<j>) 4- P 21 (s<f>) (c£ x cA 

4* s 21 sA) + P 22 (s<{)) (c 22 c 2A t s 22 s2A) 


(92) 


where the c ! . and s ! . are functions of the mutual orientation 
iD xj 

angles as follows: 


<?z 


- 1 r 

0 I 

a ,2 M' L 


A - + ^ + C - ~ (y 2 A + y ,2 B + y ,,2 C) 


] 


= i r 

21 a 1 2 M' L 


ayA + a'y'B + a ,, y ,, C 


i = 1 [ 

1 a ' 2 M* L 


y3A + y’3'B + y' '3.' *C 


4 

4 


(93) 
(cont 'd. ) 
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C22 = 


4a 


f (g 2 - a 2 ) A ■ + (3 * 2 - a 1 2 )b + (g” 2 - a ,,2 )C 

' 2 m' L J 


>22 


1 r 

a 

» 2m’ L 


2a' 2 M 


gA + a'g'B + a ’ ' g 


' *cj 


The force acting on the Earth projected on the ‘inertial axes 


{X^} is 


4 X? 

x 


= m. 


[ 


3U 


45 


3r 


(t 


->• * 1 
I . ) + ~~~T 

i rc<j> 


3uJ 5 

IT 






1 3^45 + 

r “3^~ (l <j) 




(94) 


where r, §, X and the associated unit vectors are the spherical 

T 

polar coordinates of the Earth with respect to {y\ •}. 

The appropriate partial derivatives are 

3uJ 5 


3r = -3Ga 2 r~ 4 j c 2 o ?2 o (s$) + P 2 2 (s<i>) (c 22 c2A 
+ s 22 s2A) | -3Ga ,2 r“ 4 | c 20 p 2o( s ^) 

+ P 21 (s4>) (c 21 cA + s^sA) 


+ P 22 ( s <|) ) (c 22 c2A + s 22 s2A)| 


(95) 


3U 4 5 


= Ga z r 


2 r —3 


| C 20 P 20 + 2P 22 (-C 22 S2A 
+ s 22 c2A j- + Ga ,3 r~ 3 | c 2 


; 20 p 20 


( 96 ) 


+ P 21 (-c 21 sA + s 2 icX) 


21 1 


2P22 ( ^ 22 ^^^ 
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auj 5 

3 $ 


- Ga 2 r~ 3 | C 20- P 20 ( j ) + P 22 ( p( c Z2 c ^^ + S 2 2 s2 '0| 

+ Ga ,2 r“ 3 |c2oP20^ + ^ 21 ^ (c 2 ic2A + s 2 is2X) (97) 

+ P 2 2^ ( c 2 2 c2 ^- "P S 2 2® 2 ^-} ^ 

Torque on Earth Due to Moon and Sun . Again considering only 
second degree terras this torque is 

M yl = 3GM i r ls“ 3 (° ~ B > m 0 n 0 
+ 3GM 3 r 12 f 3 (C ~ B)m^ n^ 

M y2 = 3GM iri5 - 3 (A - C)£ Q n 0 . 

(98) 

+ 3GM lt r 1 i t “ 3 (A - C)^ 

M y3 = aGM^is' 3 (B - A)£ 0 m o 
+ 3GM^r ^ i^ 3 (B - A)£^ 


CONCLUDING REMARKS 

This report presents a unified development of a physical model 
and a mathematical model of the Earth-Moon system. The Earth and 
Moon are considered to be rigid bodies. The equations of motion 
are formulated in a completely coupled fashion and the mutual 
potential of the Earth-Moon pair is incorporated in the development. 

This model is intended as a basis for a more inclusive theor- 
etical model including relativistic, non-rigid, and dissipative 
phenomena . 
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The models are being coded for use in data reduction packages 
to estimate physical parameters of the Earth-Moon system. The 
listing for two programs that have been developed to date are 
provided in Appendix C. 

Program ANEAM0 evaluates (1) a truncated form of Brown's 
lunar theory, (2) Eckhardt's theory for lunar physical librations, 
and (3) Newcomb's expressions for the rotational motion of the 
Earth. Program RIGEM numerically integrates the rotational motion 
of the Earth and Moon and the translational motion of all the 
planets. More information on these may be found in the listings. 
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APPENDIX A 


FORCES AND TORQUES BY VECTOR-DYADIC METHOD 


Reference 16 provides a derivation of the appropriate equa- 
tions. These are summarized below as applied to the forces and 
torques on the Moon due to the Earth. 


A. Force on Triaxial Moon Due to Spherical Earth 


GMM r i r 3 . GMQ 

2 2 u 


(Al) 


3GMi 


(E0 - I) - 


15 

2 


GMx 


-y 

r 


(E0 - I) 


-t 

i 


B. Force on Spherical Moon Due to Oblate Earth 


-t 

i 


F = - GM’M 


i 2 

+ 1 


5(J 3 * I r ) 33 


+2(^3 • ? r ) 3 3 


]| 


(a 2 ) 


C. Torque on Moon Due to Spherical Earth 


£ -3GM ->■ = -*■ 

T = l • I x x 

^3 r r 


■ (A3) 
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D. Torque on Moon Due to Oblate Earth 


$ = ~3GM 


■r 

x 


X 1 


5GMJa 


i 2 


[H- 


V(J 3 



I x x. 


(A4) 


+ 2 j 3 • x r (^3 • I x i r + x r • lx 3 3 ) 
2 '■+ = -t "I 

“ 5 13 ‘ 1 * 13 J 


In the above equations. 


0 =* (A + B + C)/2, 

“ "J" *t* "h “h 

1 = A^iJi + Bj 2 j 2 + C J3 D3 
E = J 1 D 1 + 3232 + 3 33 3 
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APPENDIX B 


DERIVATION OF TORQUES FROM THE MUTUAL POTENTIAL 

Consider the {y^} frame as the reference frame. The 
potential at any point of the moon (y lr y 2 , y3) is given by 
<j> {y lr y 2 , y 3 ) in its most general form. The total mutual 
potential may be written as 


J 4> (y 1 / Yzr y 3 ) dM ' 


Note that the subscript 45 on U has been omitted here. 

Now, the force on a particle of mass dM' at point (y l7 y 2 , Y 3) 
resolved along the {y^} axes is 


where 


-y -v -v -v 

£ - f yi 31 + f y 2 32 + £ ys 32 


components of this force along the axes {z^} are 


{f z } = a* 3 ' y* 

i 

_a" B" y" 


{f •} 
Y -i 


The differential torques about the {z^} axes produced by 


these forces are 


m . = Zof ■ ~ Z qf 

Zl z Z3 Z 2 

m = z q f - Z 1 f 


m = zif - Zof 
Z -3 1 Zz ^ Zl 


59 



These may be written as 
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then 





ayi r 

3 ^ dM ' - a 'J 




9 yi 


Y i 3a 


rr dM’ 


+ B ' * / <f> 


/ 


3 y 2 


K y 2 3 S 


dM' 


- 3 '/ 




3 y 2 


y 2 9 B 


rr dM* 


3y 3 


3 y 3 


/ d y 3 r 3y 3 

- *■/ * y 3 w^ m ' 


(B7) 


etc. 


Finally, 


M = a ’ * U 1 
z i a 


a * 



+ B* » 



B* 





M^ = a U 1 , , - a * * a 1 
z 2 a' ' a 

+ B Ug,, - B*' u* . 
+ y ^ii ~ y ' 1 


(B 8 ) 

(cant’d. ) . 
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and 


M = a' U 1 - a U 1 , 
z 3 a a ' 


+ B’ U* - 3 U* r . 


+ Y* u y - Y U y» 


(B8) 
(concl 'd. ) 
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APPENDIX C 


PROGRAM LISTINGS 



PROGRAM RIGEM ( INPUT, OUTPUT ,TAPE2= INPUT, TAPE3=OUTPUT) 


PROGRAM RIGEM 

THIS PROGRAM INTEGRATES THE TRANSLATIONAL MOTION OF. THE 

SUN, PLANETS, AND MOON . IT ALSO INTEGRATES THE FULLY COUPLED ROT ATI ONAL 
MOTION OF THE EARTH AND MOON USING SUBROUTINE RA19S 


? 

t 

i C VARIABLES 

1 


I 

m c 

— 1 

c X VECTER 


X VECTER OF COORDINATES 
V VECTER OF VELOCITIES 
XU-33) TRANSLATIONAL 
X (34-37 ) EARTH ROTATIONAL 
X (3 8-41 ) LUNAR ROf AT IONAL 



_ MULTI-CASE OPTION ■ _ 

7 ^ “17 ICODE=0 LAST CASE 

" ICODE= l RETURN FOR NEW CASE 

NO tE- “CURRENTLY CODED TO READ NEW BETA T ft.~ 
'RATES ONLY 


PRIME 



”000003 
_0 0 0 0 03 

000003 
_000003 
000003 
000003' 
000003' 
>00004 
000006" 
0000 L0"“ 


D I MENSION X{41),VI4L) , NORO (2 ) , OMPL < 1 ) ,D{4,4) , FI41) " 

DIMENSION T ( 3 » 3 ) » E t 3 , 3 ) ? C { 3 , 3 ) ,Pi3,3) t R ( 3 , 3 } ~> P P~( 3 DEL~( 3 ) 

_1 7 DAL (3) » DELV { 3 ) , DALV{3)', PO { 3, 3 ) , ED { 3, 3 ) ’ ] 

_ D I MEN S I ON 0MM(4) t 0MDM14) ,RT(43 

_ INTEGER OMPL ' 

_ C GMMON/PARAM/ NO R 0 , AO , A D , V JJ3E P, 0 MPL , TJ NC_, T MAX 

CCMMCN/XGUT/OMM , OMDM ,_RT 

1 1 COD£ = 0 1 — 

">1^=3.14159265358979 ~ " 

_>TR = PI/180. " ~~ 1 1 

RTD=i 80 • /PI 


GET INITIAL CONDITIONS 


000011 


11 CALL "PROB ( X t'V »T I , TF , N V , NC LASS , NS S, N I » NOR, LL » II CODE ) 


INTEGRATION 


4 CALL RA 1 9 S (X » V , T I , TF »"X L , LL~» NV,N I , NF , NS V NCLAS S , NOR, NS S) 


000024 




o' o jo o'o o o lob'oo o'o. b olo ob o‘o o o o o o o o o o'o o 

obol ; O’O O 0 : lolo’o O O O j o 0-0 Oio OOOOOOOOOO O’O o 
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M.N) ro, 'J\J Mjfv t\> ho 1 !" I— r— t— (— I ,►_» >»— it— I—.— I— H-l— OOOOOO 

o-;ui x-; h-" ON'o> O' ui ui' i -t'-w utuj Mi-ooo->i^ui4 s 4'j' 

^-,U> I 'Oi^|tjJ|-J, UIIOJ -sj O' Ul i^-'iO'O' I— I— 

111 ! [ lit. .III... . i . i 


c 

c 


OUTPUT,. SECTION 


TWR=VJDEP+TF . 

"_W I TE 1 3 , 100) 7 _7_ 

____ WR I TE (3,115) TWR 

_ WRITE (3, 101 )VJDEP,NOR,LL, NF, TF 

WRI TE (3 » l 05 ) 1 

WRITE<3»106) 

I 1 = 1 

DO 1 1=1,31.3 

_ WRITE (3,102) I I 

WRI TE(3, L03)X( I ) » X< I + i )., X ( 1+2 > 

V( I )— V( 11*100, _ __ 

V( I+1)=V( I+l)*100. 

V(I+2)=V( I+2)*100. 

WRITE(3tl04lvliLtVfjilJbV(J ±2 L 

11=11+1 

_1_ CONTI NUE 

_ DO 75 1=1,33 

75 V < I.ifVU) / 1 0 0 , 

C . . E A RT HJQRI 

IF (NOR 0(1) . EQ »0 ) ~G0 

WRI TE(3,107) " ^ 

WRITE (3,115) TWR' 

WRITE(3,108) 

WRITE (3,111) ( X ( 1+33 

WRI T EJ 3,111) { V < 1+3 3 

C _ _ 

C _ _ EULER_ PARAMETER TESTS FOR EARTH 

c ' ' 

WRTtE (3,113) _ ■" 

JJE S T 1=X ( 3 4 ) tt* 2+ X ( 35}**2+X( 36 )W2 + X (3_7)**2 

T E s7 2=X < 3 4 V* V < 34 3 +X ( 3 5 ) * V ( 35jT+ X( 36 3*V(3 6)+X( 37)*V (37) 

WRJ TE (3 >.112) _TEST1 , TEST 2 . 

C "”"77" M 0 O N "OR I ENT AtY6 N 

76~H RI TEj (3 v 10 9~) ~~ 

WRITE ( 3, 115) JTWR. 

WRITE (3, lTo ) 



) , 1 = 1 , 4 3 _ 
) , 1 = 1 , 4 ) 




000265 _ WRITE(3»111)(X( 1+37) ,1=1,4) 

_00Q2 77 WRITF{3 ,111) ( V ( 1+37) , 1 = 1,4) 

_000311 _ WRITS (3,116) 

__ _ " C “ ' 

__ __ £ _ ROUTINE TO CALCULATE EARTHS SEL ENQGRAPHIC COORDINATES 

_ _C 

000315 D( 2 , 2) =X{ 38 )»*2+Xt 39) **2-X (40 )**2-X (4 1)*»2 

000323 D(3,3)=X(38)**2-X(39)**2+X(40)**2-X(41)**2 

000331 . D(4,4)=X(38)**2-X(39)**2-X(40>**2+X(4L)**2 

000336 J D { 2 , 3 )-2 • * ( X ( 3 9 ) * X ( 40 ) +X ( 38 ) *X ( 4 1 )) 

000342 ~ D { 2 , 4 ) =2 . * ( X ( 39 ) *X ( 41 )-X ( 3 8 ) *X ( 40 ) ) 

000346 0 ( 3 » 2 ) = 2 . * ( X ( 3 9 ) * X ( 40) -X ( 3 8 ) *X ( 4 1 ) ) ~ 

000352 D ( 3 , 4 ) = 2 - * ( X ( 40 ) *X { 4 1 ) +X 1 3 8 ) *X ( 39 ) ) 

000356 D ( 4 , 2 ) = 2 . * < X ( 39 ) *X ( 4 1) +X (3 8 ) *X (40 ) ) 

000362 f)(4,3) = 2. *{ X ( 40 ) *X ( 4 1 )-X ( 38 ) *X( 39 ) ) ' " ' 

_0003 66 DAL( i) = X( 13 ) -X (10) $ DAL ( 2 ) = X ( 14 )-X ( 1 I ) $ DAL (3 )= X ( 15 )-X ( 12)” J 

000374 RRR = SORT { DAL ( 1 ) *DAL ( 1 HD AL ( 2 ) *D AL ( 2 ) +DAL ( 3 ) *DAL( 3)} 

000402 CC=DAL(1) /RRR $ CS=DAL f 2 ) / RRR $ SPH=DAL(3 )/RRR 

000406 CPh=$QRT<DAL( 1)**2+DAL(2)**2J /RRR 

000414 _____ CL = CC/CPH $ SL = CS/CPH ” “ 

000417 CALL FORCE ( X , V , T F , F ) 

000422 ________ SS= SORT { D { 2 , 2 ) **2+D ( 3 , 2 ) **2 ) ’ ’ 1"" 

000430] _ SLONG=ATAN2(D( 3,2) ,D< 2,2) ) 

000433 SLAT=ATAN2(D(4»2)»SS) ~ ' ‘ “ 

000436 _ _ SLQNG=SLONG*RTD $ SLAT=SLAT*RT D " "" 

000440 ~ WRITE (3,117) SLONG, SLAT 

_C _ ' ------ 

1" C " _]]] TEMPORARY OUTPUT 

"C ' ' ' 

_000450 WRITE" (3,9000) ChMM ( 1) , OMM ( 2 )70MM ( 3 ) » OMM ( 4 ) ’ 

000464' 9000 FORMAT ( 5X»*CHECK* f 4E20. 12 ) 

000464] - WRITE (3,9000) OMDM ( L ) , OMDM ( 2 ) ,OMDM ( 3 ) , OMDM l 4 j 

00050 0. ' "W RITE (3,9000) RT ( 1 ) , RT ( 2 ) , RT ( 3) , RT ( 4 ) 

C ‘ 

C “ ROUTINE TO CAL CUL AT E * RHYS ICAL DERATIONS 

__c ' " _ _ 

000514 T ( 1 , 1 )--CC $ T ( L , 2 )=-CS $ T(1,3)=-SPH 

000521 T ( 2 , 1 ) = SL $ T(2,2)=~CL $ T<2,3)^0. 

.000525 T ( 3 , 1 )=— CL*SPH $ T ( 3, 2>=-SL*SPH $ T(3,3)=CPH 

000 531 TEP= ( VJ DE P+TF— 24 L5020. 0)/36525. $ T£P2 = TEP*TEP_ $ TEP3=TEP2* 

0 00537 EPS=( 23.452294 -.0130125 *TEP-. 00000164 ""#TEP247 G00000503 * 

1)*DTR 

" e < 1 7 1 )* 1 . " "$ E { 1 ,'2 ) = 67” $ E (1 , 3) =0* 



000546 


CL.rO 
LU CL 

>-;lu 


000551 E £ 2 , 1 ) =0 . $ E £ 2 , 2 I = COS(EPS) _ $ E ( 2 , 3 ) =- SIN(EPS) 

000556* E(3,U=0. $ E<3,2)= SIN(EPS) $ E£3,3)= COStEPS) ..... . 

000564 011,11^(2,2) $ C(1»2)=D(2,3) $ C ( 1 , 3 )= 0 ( 2 , 4 ) . 

000570 C(2 t l) = P(3,2) $ C ( 2* 2 ) = D £ 3,3) .. $ C( 2, 3)=D( 3 t 4 ) . ... 

000575 C<3,1)=0(4,2) $ C ( 3 , 2 ) = D ( 4 T 3 ) __ $ . C £ 3, 3 )=D ( 4,.* ) . . 

000601 TT=VJDEP+TF — 2433282. 5 _ 

000604 XKAP=. 063107 *TT*DTR/3600 . ... . . _ 

000610 OMEG= .063107 *TT*DTR/3600. . . 

000612 XNU=. 0548757 *TT*DTR/3600. ■ ' . . - 

0006 14 ' ~ PI 1,1)=- SJN(XKAP)* SHN£OMEG)+ COStXKAP)* COS£OMEG>* CCS(XNU) 

000 6 33^ " P ( 1 » 2 ) =— CDS ( X KAP ) * S IN( OMEG J— SIN( XKAP)* COS (OMEG ),*„ COS(XNU) 

000652 >( 1 ,3)=-" COS(OMEG)* SIN(XNU) 

000660^ P £ 2 , 1 ) = SIN(XKAP)* COS (OMEG )+_CQS( XKAP >* S IN < OMEG >*_CCS( XNUJ. 

000677 P £ 2 , 2 ) = COSt XKAP) * COS(OMEG)- SUHXKAP)* S IN tOMEG ) * COSl XNU) 

000715 P i 2 » 3 )=— SINtOMEGJ* SINtXNU) " 

0007 22 P £ 3 ? 1 ) = COS ( XKAP ) * SIN(XNU) 

*0007 30 P ( 3 j 2 )=— S I N £ X KAP )* SIN(XNU) 

000735 P ( 3 » 3 ) = COS ( XN U) 

000740 DO 50_.K= 1,3 

000741 " DO 51 J= 1 , 3 

000742 R ( K t J ) = 0 • 

000745 DO 52 1=1,3 

J3 0 0 7 4 6 H_DO 53 L = l, 3 1" ___ ■ 

000747 _ 53 P (K, J)=R£ K, J )+E ( I ,K )*P U »L)*T( j,U 

000764 52 CONTINUE 

000766 Si’ CONTINUE' ' _ _ ... 

000770' 50 CONTINUE . ... 


C 

c 


MATR I X PP BECOMES THE PRODU CT E(TR)*P *T( 


TR ) *C (TRJTR=TRANS 



C 

000772 DO 70 1=1,3 

“000774 DO 71 J = l,3 „ _ 

*000775* PP(I,J) = 0. 

"001000_ DO 72 L=1 ,3 

*_001001 *7 2 P P £ I , J ) = P P ( I , j ) +R j I , L )*CC J,L) 


001014 7 1 CONTINUE. 


001016 70_C0NT 1NUE __ 

001020 APHI = ATAN2(-PP£3, 13 ,-PP£ 3, 2) ) 

001026 ST= SQPT( PPC3, 1)**2+PP{ 3,2)**2)_ 

001033J ATH= ATAN2(ST,PP(3 t 3) ) 

0010 36 AP_SI=„ATAN2(-PP(1,3) , PP(2,3) ) _ 

00104 2 A PH 1= APH I *RTD t ATH=ATH*RTD $ AP S I=APS I* RTD 

001045 APH I = AMOD£ APHI ,360. *) 




001050 

001053 

"001056 

001060 

^001072 

001106 

001122 
001124 
_0 0 1 1 2 7 
001131" 
001L34 
001137 
001151 
001153 
j)01 155 
00 1161 
001166 
001172 


IF (APHI.LT.O.) APH I = APH I + 360. 

APS I = AMOD( APSI ,360. ) 

IF (APSI.LT.O.) APSI=APS 1+360. 

WRITE (3,120) APH I » ATH » APS l 

A0MEG=259 .1832 75-0. 0529539222* ( 36525. *TEP ) +0.000 1557* I 3. 6525 
. **2) *TEP2+ 0.0000 0005* (3. 6 525**3 }*TFP3 

AM00N-270.434358+13 . 1 763 96 5268* { 36525 .*TEP )-0 .000085* (3. 6525 
**2) *TEP2+. 0000000 39* (3. 652 5**3 )*TEP3 
AI=5549. 3/3600. 

AQMEG=AMOD( AOMEG, 360 . ) 

IF (AQMFG.LT. 0. ) AOMEG= AOMEG +360. ' “ " 

AMGCN=AMDD(AMQQN,360.) 

IF (AMOON.LT. 0. ) AMOON=AMCON+360. " " ' " ^ 

WRITE (3,121) AMOON,AI, AOMEG Z. 

RHO= A TH— A I __ 

S IG=APS I-AOMEG ' ' ’ 

TAU=APHI-180.“AM00N+APSI 

IF (ABS(TAU) .GT.10. ) TAU=TAU-360. 

WRITE (3,118) 

WRITE (3,119) RHO ,SIG»TAU 


001204 
001212 " 
001221 " 
001226" 
00 1236 
_0 0 1 2 4 0 
001243 
001255 
001260 
001261 
"001263 
J) 0 1 2 63 
001263"' 


001263 

001263" 

J301263J 

‘001263 


100 

ToV 


102 

103 

“lO£ 

105 

"106" 


EULER PARAMETER TESTA FOR MOON 

"TEST3=X(38)**2+X(39)**2+X(40)**2+X(4L)**2 " 

TFST4=X{ 38)*V(38)+X(39)*V( 39)+X(40) *V (40 ) + X( 4 1 )* V (4 !)”"_ 

CONl^-( V(38)*V(38 ) +V (39 j *V (39 ) +V (40 ) *V ( 40 ) +V ( 4 1 ) * V( 41 ) )_ 

CON 2=X( 3 8 ) *F ( 38) + X( 39) *F ( 39 ) +X <_40 ) *F ( 40 ) + X ( 4 1J *F ( 41 j 

DC0N=CDNi-C0N2 

WRITE (3,113) _ ______ " 

WRITE (3, 112) TEST3 , TEST4, DCON 
IF (TF.GE.TMAX) GO TO 2 
TI=TF 
TF=TF+TINC 

GO TO 4 "" " 

F0HMAT(*1*,40X »*TRANSLAT I ONAL~‘MOT ION**// ) "" 

FORMAT (* * , * E P OC H=* , E 1 9 . 1 2 * 5 X , *ORDE R« * , I 2 , 2X” *LL=* , 12~,2X ,*NF=* 
i»3X 

.»I5,*DAYS PAST EPOC H=* ,£ 15 .7 ,7/) 

FORMAT!* *,*PLANET*» 14) 

FORMAT! * 

FORM ATJ * 

“forma't ' ' 

"FORMAT" 



=k 


tu 


JHL 


*,3(E19. 12, 2X) ) 
*,3( E19. 12, 2X)//) 


(9X,*X( AU) *»14X,*Y( AU)*»14X,*Z(AU)*,/j 

T 9 X , *XD l AU/ LOOD j * , 9X ,*YD ( AU/ 10 OD ') * ,9XV*ZD< AU/lOOD)* * / T 


001263 



001263 107 F0RMAT(*i*,40X,*EARTH ROTATIONAL MOT I CN* , // ) 

001263" 108 FORMAT (20X,* EULER PARAMETERS AND RATES, BETA PRIME*,/) ... 

001263 "l09 FORMAT(*l*,40X,*MOGN ROTATIONAL MOTION*,//) 

001263 110 FQRMAT(20X»*EULER PARAMETERS AND RATES, BETA TR.PR IME*, /) . _ ... 

001263 111 FORMAT (9X,4EL9.12) ' 

001263 112 FORMAT (1 OX , 3 ( E 19 . 1 2 , 3 X ) » / / ) . 

001263 113 FORMAT { / ,20X,*FULER PARAMETER TESTS*) • 

001263 115 FORMAT <* * ,40X, *JUL IAN DA T6=* ,E 19. 12 , //) . . . . 

001263 1 16 FORMAT { / / , 20X ,* E AR THS S EL ENOGR APH I C COORDINATES*,/) ' 

001263 117 FORMAT (* * , 9X , * LONG .= * , F 19 . 12 , 5X , *L AT „ = * , E 19 . 1 2 , // ) 

001263 118 FORMAT! 35X,*LUNAR PHYSICAL L I BRAT I ONS *, / ) 

001263* 119 FORMAT ( 2 X , *RHO=* , E 19 . 12 , 5X , * S IG=* , E 19. 1 2 , 5X , *TAU=* , E 19 . 12 , // ) 

001263 ” 120 FORMAT (35X,*EULER ANGLES*, //, 5X , *PHI =*, E 20- 12 ,2X ,*THETA=* tE 20_ sl 12 

1 ,2X,*PSI~*,E20.12,//) ■ 

001263 1 21 FORMAT ( 35X, *FUNDAMENTAL ANGLE S* ,//, 5X,*M00N=* ,E20. 12 ,2X, * I=* 

1 ,E20. 12, 2X, * OMEGA*, E20.JL2,//,) „ ■ 

001263 122 FORMAT (15) _ 

0*01263 2_R E A D (2,122) LCODE ■ 

0012 71 I F ( I CODE . EQ. 1 ) GO TO 77 _________ 

001273 STOP _ ' 


001275 END 





NSF = . FALSE. 

NCL=NCLASS.EQ. 1 

NPQ=NCLASS.LT.2 

SR= 1.5 

IF(NV.EQ.l) SR= 1 . 2 

NES=LL. LT.O 

TDI F=TF-T I . 

DIR=TDTF/ABS(TDIF) . 

IF(NES) XL = ABS ( X L) *D IR 
NCLASSM ABS(NCLASS) 

LA=KC2*KD2-l 

'00 14 N=2,KF _ 

LA=LA+1 

~H(N)-=HH(LA) 

W(N-l) = ONE/FLOAT ( N+N * *2* l N C L A S Sr. ’ 

U ( N- 1 )=N+ 1 

JDO 22 K = 1 ,NV 

I F { NCL J V ( K l-Z.ERp 

DO 22 L=1,KE 

BT(L,K)=ZERO _ . . 

R ( L »K) = ZERO 

Wl = CNE/FLOAT (NCL ASS*) _ 

DO 9 39 J = 1,KD_‘ 

M=MC ( J ) _ 

JD= J + 1 “ "1 

CO 939 L=JD,KF 

XI ( P)=FLOAT (NX I (M) )*W(_ JJ/WCLI 

’M-M + l _ __ 

JI( 1 }=-H(2)*W{ 1 ) 

0(3.) -H ( 2 ) / W ( 2 ) _ 

R { 1 > = ONE / ( H < 3 )— H ( 2 >j 

LA= 1 

LC= 1 

''do 73 K= 3 ,KE_ 

_L B= LA ‘ 

LA=LC + 1 

LC=NW(K+l'f _ 

JD=LC— LA 

' C(LA)=-H(K)*C(LB) 

"C(LC)=(C (LA-1)/W( JD)-H( K) )*W( JD+: 

j) ( LA)_=H ( 2 ) *D ( L8 ) *W{ K-l ) / WJ K) 

_D ( L C }= ( D ( LA-1 ) *W (_K- 1 ) -HH ( K ) )/W (K ) 
R (’ LA ) =ONE/ ( H ( K+l )~H 12) ") 


000240 R(LC)=0NE/(H(K*1 )-H(K ) ) 

000246 IF(K,EQ.3) GO TO 73 

000250 _ 00 72 L~4,K 

000251 "" LD=tAtL-3 

j) 002 54 LE= L8+L— 4 .. _ 

000256* JDM=LD~L A 

"000260 C(LC)=W< JDM+1)*C( LE ) /W I J DM )-H( K) *C ( LE U) . 

000270 _ D<LC)=(D( LE)+H(L-l) *D(LE+1) ) *W I K-l > /W ( K ) 

'000301 _"_72 " R (LC)=ON£/( HiK + l)-H(L-in 
"0003 11 _ 73 CONTINUE 

000314 ' SS= IQ.**(-LL) ... 

000320 _ ~NL = N I +30 

_C “SET IN A REASONABLE ESTIMATE TO T BASED ON EXPERIENCE, (SAME SIGN AS TF-TII 

"000323" _ I F ( .NGT«NE$) TP= (( FLO AT( NOR ) / 11* ) *0. 5**{ 0. 4*FLCAT ( LL ) ) ) *DI R 

000337 _ __~"IF{NES) TP=XL " 

' 000342 “ ' IF I TP/TDIF,GT.0.5 ) TP=0. 5*TDIF .. . . . _ 

JD0034 7 NF = C 

000350 ~NC0UNT=0 _ _ ... . _ 

"000351' 4000" N S=0 _ ... . . . , 

"000352 _ “ -TM*TI " . . . 

000353 " " SM»L.E4 . ' 

000355 CALL FORC E ( X t V * TM ,F U ,. . __ 

"000357 " ' NF= NF + 1 f _ _ 

C LOOP 58 FINDS THE" BETA- V ALU ESjFRQM “THE CORRECTED B-VALUES, USING p-C EFF 

000361 722 DO 58 K=1,NV 

000366' BE( KE,K)=B(KE,K)/W( KEJ 

000374 DO 58 J=1,KD 

"000375 ~~"jD=J*i 

"000377' BE( J«K)/W( J) 

000404 DO 58 L=JD,KE 

000406 "N-NWIU+J ’ 

000411 58" BE(J»K)=BE(J»K1+0(N)*B(L»K) 


00 0431 T=T P _ _ ... 

000432 "" TVAL= ABS ( T ) "' ' ' 

000434 T2= T^^NCL ASS “ ' _ _ _ _ _ 

C LOOP 175 IS THE "ITERATION LO“OP_ WITH_NL=NI PASSES AFTER THE FIR ST SEQ E NCE 

000440 " ““ ' DO 175 M-l f NL* " 

000442 J2= .TRUE , ' 

000443 DO 174 J*'2yKF I 

000445 1”"J0=J-1 " " ________ 

0 00447 L A= N W ( J DJ " 

000451 J DM- J -2 


0i00453 S=HU ) 



000454 Q*s**(NCLASS-U 

000462 IF(NPQ) GO TO 5100* 

000464' DO 1300 K=1,NV 

000465 ' ~ RES=6(KE ,K) _ ' 

J300470 TEMP=RFS*U( KE ) 

000473 DO 7340 L=1*KD 

000474 _ J R= KG— L 

000476* _ RES = B( JR,K)+S*RES 

000503 ^ 73 4 oj* TE M P=B ( JR , K ) *U ( JR ) + S*T£M P 

000512 Y{K) = X(K)+Q*(T*V{K)+T2*S*(Fi{Kl*Wl + S*RES"V) 

10005 3L 130 0 ' Z {K)=V<K)+S*T*(F1(K) »S*TEMP J ™ 

000543' _ * __GO'lO 5200 1 

0005 43 "5 100 DO 1400 K = l f NV " HI ’ ~ 

000545 __ _ RE$ = B{ KE» K > 

000550 DQ 2340 L = 1 , KD __ 

000552 JR=K?-L ‘ J • 

_0 005 54 2340 R E S= B { JR , K ) +S* R E S ’ H 

0 00 5 64 1400 Y { K ) = X { K j +Q * { T * V ( K )_+T2* S *j F 1 < K )*_W 1+ S* RES ) ) 

000604 52 00 CONTINUE 

*000o04__ CALL FORCE! Y»Z,TM+S*T»fJJ _H _ H 

1000612 'NF=NF+1 

000614 IFU2) GO TO *702 


DC 471 K= 1 1 NV 

"TEMP-=BE(JD »K) 

R E S= ( F J {K ) - F 1J K ))/ S’ 

N=L A _ 

DO 134 L = 1 » J D M* 1_ 

N=N+ 1 

134 RES= ( RES-BE ( L » K 1 

BE{ JD tK) = RES 

TEMP=RES-TEMP 1 

B(JD.K) = B(JD,K)+fEMP*W(JDJ 

N=L A 

" DO 471 L= l f JDM _ 1 ' 

1 N«N*1 _ _*__ 

471 B ( L , K ) = B t L » K ) +C(N )*7EMP 

GO TO 174 

70 2 J2=. FALSE. " __ 

D O 271 K- 1 1 NV ~ 

TEMP = BE( 1 »K ) ~ 

R ES= f FJ (K i- FH K) ) /S 

BE ( i * KJ-RE S __ 

2 7 1 B (I , K ) = B { 1 , K ) + C R ES-T EM P 3 *W I 



000724 174 CONTINUE 

000727 21 IFf N.LT.NI1 GO TO 175 

000 7 31 HSUM=0. 

‘000732 __ VAL=TVAL**(-KE) 

” OOP 737 00 635 K=l t NV _ 

000741 635 HSUN'=HSUM + B(KE»K}#*2 

‘000751 2 HS UN=VA L* SORT ( HSUM) 

000754 2 ” IF C NSF ) GO TO 175 

000761 “ ‘ I F ( AB$( <HSUM~SM)/HSUM).LT*O.Ol) GO TO 176 

'000766 ” " SM= HSUM 

”000766 1 7J5 _ CONTINUE " _ 

C' THIS NEXT PART FINDS THE PROPER STARTING VALUE FOR T 

’000771” 176 I F { NSF ) GO TO 180 

000773 " ■' IF(.NCT.NES) TP= ( SS/HSUM ) **PW*DI R 

001002; IF(NES) T P= XL __ _ 

”001005 I F< NESI GO TO 170 7 __ 

200L006; 2. IFCTP/T.GT.ONE) GO TO__170L 

001013 8 FORM A T ( 2X1 2 1 2E 18 « 10 ) _____ _ _1 " 2 

001013 “ ' TP = 0 . 8 *TP " “ 

001014 'NCCUNT=NCOUNT+L 

"001016” I F t NCOUNT »GT* 10 ) RETURN 2 ‘ 2 21 1 " 2 " 

001021” ” ‘ PRINT 8 » KD » T t TP " " 

001033 GO TO 4000 

001037 _____ 17C PRINT 8,KD,T,TP" ' ” " _ __ 

0010 51 N SF = * TRUE . 

C "FIND POSITION < AND VELOCITY FOR CLASS II “AND ITS) AT’THE END OF THE SEQUENCE* 

" 001052 ' 22180 _ DO 35 K=1»NV 2 _ __ ' 

001057 ’ 2 RES=B(KE»K) _ ''2" 

' 00L06 2 " ' " DO 34 L=1»KD 22 2_2 12 

001064 34 R£S=RES+8 ( L» K ) ~ _ 

001072 __ _ X(K)*X<K)+Vm*T + T2*(Fl<'lO*Wl+ftES)2_ _2 

”001105 ‘ I F ( NCL ) GO TO 35 ' “ " 

001107 RES=B(KE,Ki*U(KE) 2 __ ” 

001113 ‘ ‘ DO 33 L=1»KD 

00111 5" ' 33 RES=RES+B<L ? K)*U(L)"~ " 

001126 _ 2 V(K) = VIK)+T*IF1(K>' +RESJ __ 

‘001133 ’ 35 CONTINUE 2 2 '"22222 _ 2 _ ' " 

001136 ’ TM= TM+T 2 2 2” ’ " 

” 001140 ” 2 ns=ns+i ‘ “22 ' 1 ’ 22 

001141 74 l F { NPER ) 'RETURN ' ”'“ 

C CONTROL ON SIZE OF NEXT SEQUENCE AN‘D‘”RETURN'“WHEN TF'"lS REACHED 
001144 CALL FORCECXy VvTMy'Fl) 

001146 NF*NF+1 



0 01150 IFINES1 GO TO 341 

001155 TP= <( SS/HSUM)**PW )*OIR 

001163 _ _ IF { TP/T.GT.SR) TP=SR*T 

001167 341_IF(NES) TP=XL 

001172 IF{DIR?(TM+TPKLT.DIR*TF_ri.£-10) GO 

_001200 TP-= TF-TM 

001202 NPER = „TRUE* " ' '' __ 

C PREDICT B-VALUES" FOR NEXT SEOUENCitl 

001203 77 0 =TP/T _ 

001205 DO 39 K=l f NV 

001207 _‘res=cne ' ‘ "~7‘~ H _ 

001211 , DO 39 J-l , KE 

0012 12_ IF ( NS .GT * i ) BT ( J » K ) = B ( J , K J-BTI j t KJ._ 

001222 IF( J.EQ.KE) GO _T0 740 __ _____ 

001224 M=MC(J> _ 

0012 26 J D = J+ 1 ‘ ““ I ' "Z II 

_ooi23o’ Ido ao l=jo»ke "~ H. 

00123i_ B( J»K)-B( J t K)+XI (M)*BXL f JiI __I 

00 12 41 40 M-M+ 1 

0012 45 740 RES = RES*Q _I 

001247 T EMP-RE S*B1 J y K ) ' II 

_Q01 2 53 B(J t K)=TEMP» BT (J,K) I__ 

_00i2£.0 39 B TI J ? Kj = TEMP 

_0 6 1266 NL=NI 


001 2 67 G 0 TO 722 

001270 E ND 




McJ 


-w§- 

fh ra 
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" SUBROUTINE PROB ( X » V, T I , TF , N V, NCLA S S , NS S , NI ,NOR , LL , 1 1 CODE) 

c ’ ™"L“ 

C ~ " THIS SUBROUTINE PROVIDES CONSTANTS AND 'INITIAL ” 

C_ CONDITIONS FOR RIGEM ” 

C NORO ( I ) = 0 NEGLECT r ARTH ' ROTA T I GN~ ’ * ___ 

C “ ' =1 COMPUTE EARTH ROTATION ” “ 

C NORO ( 2) -0 NEGLECT MOON ROTATION 

C ' ' =1 COMPUTE MOON ROTATION 

C ' ' ■ - - -j”-— 

C~ OMPLI 11=1 OMITS MERCURY* SATURN , UR ANOUS, NEPTUNE , PLUTO 

c 0MPL{I)=6" 'OMITS NO PLANEtS 

c" ' 1" 1 ’Z_ 1 

DIMENSION XMASS(ll) » X 1 1 ) , V < 1 ) » XMAU 1) , NORO ( 2 ) «0MPL(1| ' 

DIMENSION XSI41) , VS(4i> " 

INTEGER OMPL 

COMMON/XMAS/XMASS 

"COMMON/ 1 NERT/BEM, GAM,'ALPM, BEE, GAE* ALPEt El 1* E I 2 * E 1 3 » AM 1 1 , AM 12 * A_M I_3 

COMMON/ PARA M/NORO»AO,AD*VJDEP»OM PL ,TINC,TMAX 

DATA XMA /.L, , 5983000. T 40 8 522. t 332945. 56 19 25 A A *27068807* 1301» 

13098700. » 1047. 3908, 3499. 2, 22930., 19260. , i.812'000./" 

IF I I ICOCE .NE .0) GO TO 12 ' "" “ 

READ 12,120) NOROIl ) ,NORO( 2) ,0MPLI1) 

READ (2,121) VJDEP, TINC, TMAX 

' REAP (2,123) Nl , NOR, LL 

PI = 3. 14159265358979 

D TR=P 1/180. $ RTD-180./p‘l 

BEM=. 00063 $' G AM=. 0002 $ ALP M=. 00043 

ALPE=. 00322 ' " " " 

3EE=. 00327 ' ' 

“GAE=. 000054 

A 0=100. 0755 42 *DTR '' ~ 

A 0=360.9 8 564 7348 *DTR 
EI1=5. $ E 12=5, $ E 13=5 . 

AMI 1=5."$ AMI 2=5, $AMI3«57" 

C THIS LOOP PROVIDES PLANETARY I.C. PER TABLE 10 , P.274, CEL. MECH. 

C JOURNAL V5, N0.3 " 

XK2=( .01720209895)^*2 

DO 1 1*1,11 ' “ ' 

L_XNASS( i )=XK2/XMA II) 

C THIS LOOP 'PROVIDES"' PLANETARY l.C. PER" TABLE" I6,p".274 

KK= 1 ~ ' 

WRITE {3 , 102) “ ' 


000124 



WR I TE ( 3 , l 03 ) 

000130 



UR 1 TE(3. 104) 

000134 

000140 

— 

— 

WRITS (3* 105) _ . 

WRITE<3, 106) 

000144 



00 2 1=1*31.3 

000151 



READ (2. 107) X( N . X ( I +1 » . X ( I +2 ) 

000176 



xsm=xm $ xsii+n=x( i + n $ xst i+ 2 )=xn+ 2 ) 

000210 



READ 12. 1071V f I ) .VU+IKV (1+2) 

000236 

000271 


— 

WRITE13.100) KK,XUi,X(IH)*X(It2) 
WRITE (3, 100) KK • VM) • V( I+l). VU+2) 

000324 



vm = vm/ioo. 

000332 



v( i + i j=v< r + n/ioo. 

000333 



V( I+2)=V< l+2)/100. 

000335 



VSIIHIII $ VS( I+1WV1 I + lt $ VS( I+2)=V( 1+2) 

000343 



KK=KK+1 

000345 


2 

CONTINUE 


C 


"THIS LOOP PROVIDES INILIAL EULER ANGLES AND RATES FOR EARTH MOON 

000347 



WRI TE(3,108) 

000352 


108 

FORK AT (/ »30X, DEARTH INITIAL EULER PARAMETERS AND RATES*,/) 

000352 



DO 3 L=34,38,4. 

000357 



READ( 2, 101) X(L) * XI L + 1 ) » X ( L+2 ) ,X < L + 3 1 

000412 



XS(U-XU) $ XS(L+l)=X(L+i) $ XS <L+2)=X(L+2) $ XS (L + 3 ) = X(L +3) 

000426 



READ ( 2, 10 1) VI L) , V( L+ 1 ) » VI L + 2) * V( L + 3) 

000462 



vs < l ) =v { i) $ vsa+n=va+n $ vsa+2)=v(L+2) $ vsa+3)=v(L+3) 

000476 



WRI IF (3, 109) X(L) fXlL+1) , X ( L + 2 > , X < L+3 ) 

000532 



WRITE (3, 109} V(L) ?V(L+l),V(L+2)»V{L+3) 

000572 



I F ( L.EQ.38) GO TO 3 

000600 



WRITE{3»110) 

000603 


3 

CONTINUE 

000611 



T 1-0 . 

00G611 



TF - T I NC 

000612 



NCL ASS=+2 

' 000613 



N V- 4 L 

000614 



NSS=0 

000615 



WRITE (3,125) 

000621 



WRITE (3,126) 

000625 



DO 11 1=1,11 

000632 


11 

WRI TE (3, 127) I , XMASS( I ) 

000647 



WRITE (3,128) ALPE,QEE,GAE 

000660 



WRITE (3,129) ALPR, 8ER, GAM 

000672 



WRITE (3,130) VJOEP,TINC,TMAX 

000704 



WRITE { 3 , 131) NORO(1)»NORO(2) »0MPU1) 

000716 



WRITE (3,132) N I , NOR, LL 


__c _ _____ _ 

"" _I c ™'~MULtI-CASE SEGMENT” _■ ' 

C 

000 730 IF ( IIC0DE.EQ.0) GO T0_14 _ _ 

_0007 35 12 DO 13 1=1,41 " . 

000737 X( I ) — X S I I ) _ 

'000741 13_ V{ I )=VS< I ) ” ' ‘ “* 

'000745 READ (2,101) V ( 38 ) , V ( 39 } , V ( 40 ) , V ( 41 ) ” ___ _ 

000 775 WRITE (3,133) _ _ 

001001 ' WRITE (3,110) ‘ ” ” 

_0 0 1 0 0 5 (3,109) V(38) ,V(39) #V(40),V(41) _ 1 

001041 “ T 1=0 • 

001045 ' TF=T I NC ’ ' 

001046 "lOO FORMAT { 3 X , 1 6 , 3E25 . L 4 ) _ “ 

‘001046J 101 FORMAT (4E20.12) ’ 

JJ01046 102 FORMAT < * 1* , 40X, * I N I TI AL CONDITIONS AND PARAMETERS*,//) 

001046 103 FQRMAT~(3GX»* INITIAL POSITION AND VELOCITY OF PLANETS*, /) 

001046 l04 'F0RMAT“(30X, PREFERRED T'O MEAN EQUATOR AND ” EOU UMOXjOF,' J?5p*_0*, //) 

_001046 105 FORMAT ( 2X , *PL ANETS* , 1 6X » *X* » 2 5X , *Y* , 25X , *Z* j _ 

001046 106 FORMAT' < 2 5X »*XD* » 25X , *YD* , 2 5 X ,*ZD* , /) ' ' 

001046 107 FORMAT (3E25.14) 

"001046 109 FORMAT ( 20X , 4E20 .12 , / ) “ ~ 


FORMAT ( 3 OX , *MOON INITIAL EULER PARAMETERS AND RATES*,/) 

"FORMAT '(3 15) ’ ” " " 

FORMAT { 3E20. 12) 

FORMAT (3 15) ” 

FORMAT <*1*,50X,*PARAMETERS*»///) 

FORMAT (40X ,*PLANET ARY GR AVI tY PARAMETERS*, / ) 

FORMAT (45X, 12 ,2X,E 19.12) 

FORMAT ( / / , 4GX p *E ARTH INERTIA RATIOS*, / ,1 5X, *ALPHA=* , E 19 . 
.* 6ETA=*,E 19. 12,3X,*GAMMA=*,E 19.12, / ) 

FORMAT ( 40X , *MOQN INERTIA RATIOS*,/, 15X ,* ALPHA=* , F19 . 12, i 
. *BETA=*,E19.12,3X,*GAMMA = *,E19..L2p/) 

FORMAT ( 10X,*EP0CH=*,E19.12,5X,*TINC(DAYS)=*,E19.12,5X ? _ 

. * TM A X ( D AY S )-=* , El 9 .12,//) ' " ' 

FORMAT (35X,*OPTIONS*,/,2'OX',*NQ'RO( 1 j'=*7l 2", 2X , *NQROl[ 2 ) =*, 
.12, 2X,*0MPL( 1)=*, 12,//) 

FORMAT (20X,*INTEGRAf ION' PARAMETERS*, /,20X,*NI=*, I3,2X,~*i 

I3,2X,*LL=*, 13,//) " " 

FORMAT ( * 1*, 50X, *NE W CASE*')*” 

iicgde=l‘” ' 


— 

&TW - 

f— 1 _ . 


B £ 

©-K- 


S 


00 1050 RETURN 


001030 END" 



SUBROUTINE FORCE ( X » V , T M, F ) 

C_ THIS SUBROUTINE PROVIDED N BODY GRAVITATIONAL “J 

C NOTATION 

C XMASS ( 1 ) = SUN _ XMASS(6)=MARS 

C ^ XMASS ( 2 5=MERCURY X MA SS ( 7 ) = JUP I TER 

__C XMASSI 3) -VENUS XMA SS ( 8 ) = SATURN 

_C XMASS (45-EARTH X MA SS ( 9 ) =UP ANQ US 

C _ XMASS ( 5) =MOON XMASS ( 10 ) =NEDTUUE 

_C" _ XMASSt 1U=PLUT0 

00000 7 ' ~ ~ " DIMENSION X(L ), V ( I) , F ( LI , XMASS < L L ), Y ( 66 ) 

000007 DIMENSION OM ( 4 ) , B ( 4, 4 ) , 8D < 4 , 4 ) ,OMD(45 , FO ( 4 ) , CMM( 4 ) »BB(4.4) 

IB8D ( 4 ,4} ,DALV(4),D(4,4), ‘ OMDM { 4) , DD( 4 , 4 ) , FT ( 4 J » RD < 4) , 

INQRO ( 2 ) ,PT<4> ,DEL<4) ,DAL<4),C(3,3> ,GMPL( U 

000007 INTEGER OMPL 

000007 ~ ^.COMMON/XMAS/XMASS 

0000 07 ’ COMMON/ 1 NERT/BEM, GAM, ALPM, BEE, GAE, ALP E, F 1 1» El 2 ,EI 3, AM 1 1 , AMI 2j_AMI 3 

000007 CGMMGN/PARAM/NGRO, AO»AD, VJCEP»QMPL»TI NC»TMAX 

000007 COMMCN/X OUT /OMMt OMDM »RT 

000007 ' R ( X I ,XJ,YI,YJ,ZI , ZJ) = SORT ( ( X J-X I) **2 + ( Y J~Y I) **2+\ 1 J-Z I ) **2 ) ’ 

000035 " C- L ♦ " 

"000037 “DO 10 1 = 1,33' " * 

_0 00040 Y ( I ) =X ( I J ' ’ 

000042 10" CONTINUE ' I 

C 1 1 IS PLANET COUNTER " " "" "" 

000044" 'DO 60 1=1,41 

000045 "60 F ( I ) = 0. “ 

000050 IF ( OMPL ( 1 ) « EQ. 1 ) 61,62 ' " ' 

000055 6 1 NN 1= 7 $ NN2=19 $ H=3 

000060 ' GO TO 63 

000061 62 NN 1 = 4 $ NN2=31 $11=2 " ' 

000064 ^ 63 DO 2 I=NN1,NN2,3 ' " 

000066 " RR=R(YU> ,Y( I) ,Y(2I ,YU+1) ,Y(3),YU+2)) 

000101_ F ( I )= -G*(XMASS(1)+XMASS( II )) *YU! / { RR*RR*RR ) 

0001 11 JJ=2 ‘ ' ' “ 

000112 I F ( OMPL ('l } • EQ • 1 ) jj=3 

000115 G1=0. “ ~ 

000 116 " G2= 0 . ' " “ 

000117 "DO 3 J=NN 1 , NN2 , 3 

000 L 2 1 "IF ( II -EQ. JJ )G0 TO 33 ' “ ‘ 

000123 G1=G*XMA$S( JJ )*( Y( J )-YU J ) /R ( Y{ I ) , Y { J ) , Y( I + 1 ) »Y( J + l > , Y{ 1+2 ) » 

NY( j+2)')**3+Gi ' "" “ _ "" ” “ ' “■’* 

000146 G2=G* XMASS ( JJ ) *Y( J ) /R( YU) ,~Y( J ), Y ( 2 U Y( J + 1 ) ± Y (' 3 J , YJ J+ 2~) ) ** 3 +G2"~ 


I i ■ j i i i * I ill ' 1 | 1 r I ! i ] ■ j i ' > I i r t t [ i i 

O O iO O O O'O O O' I O'O O O of o'o O O'O O'olo’O'O-O’O' jo O O OlQ O o o o o o 

O o lo O'O O O'O o j 0!0 ooo OOOOOO O'O oooo. ooooooooooo 

O O OOOO OiO o I 0,0 OOO jo o O'O 0:0 o 30000 jooooooooooo 

;.p- ^ -P- j ! Oj U) oj W'U> uj w ui i\j ro r\j iMWiyfyiMMryi-t- 

UltOl l\J K) • H- O' ' Oj-j — J Ul *— 1 — ' t— .t“ f—'Oi— J "J -P*. ! fSJ h- t— i-“ 1 H H O ~J *J O' 

Vnjt - 1 | W j O ■ — ■ sj’ j Wj-vj U1 W H ' Wif - . -^| ■ O' ( VJ1 Pv if-* VJl 04 r— if— 1 VJ1 -£• Cj O —4 W h* *,-*J 


33 JJ=J.J + 1 

3 CONTINUE ■ „ 

J"" F( II = F< l ) +GL-G2 

F(IH)= — G* ( XMASS ( Li +X MASS ( II) J^YC I + i ) / { RR*RR*RR J , . . _j 

J J= 2 , 

IFtCMPLli >.EQ.l) J J=3 _ ... 

G 1=0. __ _ 

G2=0 . _ ... 

DO 4 J=NN1,NN2,3 

I F ( 1 1 « EQ . J J ) GO TO 44 

_G1 = G1+G*XMASS( JJ) *_( Yl J+1J-Y( 1 + 1} )/R(YU )» Y(J),Y( 1 + 1 ) ,V { J +1 ) , Y ( I~+ 2 ) 

B,Y(J+2})**3 ' __ 

__ G2=G2+G#XMASS { J J ) *Y ( J+ LJ / R < Y (1 ) , Y { j ) , Y C 2) , Y ( J + 1) , Y ( 3 j , Y (J+_2) j **3 

44. JJ=JJ + 1 _ _ _ 

4_ C0N TI NUE " '777771 7 ~ 

F ( I + 1 ) = F I I + 1 } +G 1-G2 7 ~ 

F ( I +2 )= ' - G*<XMASS ( 1 ) +X MAS S ( Tin JY ('l+2j'/iRR_^RR^R R ) _ 

J J= 2 ’ “ ~ J ~ " 7 __1 

I F ( CMPL (1 ) . EQ. i ) “~JJ=3 1 

G 1= 0 . _ __ _ 

G2=0 . " 

DO 5 j*NNL,NN2,3 _ _ ■ 

I F ( 1 1 .CQ . J J ) GO TO 55 "7"77L _ " "" J' _ ’ _ ” I' 

__ G 1=G 1 +G* X MA S S (' J J )"* ( YJ J +2T-Y" ( ] 1 + 2 ij/RI Y~( 1 ) ,Y(_J f» Y ( I + i f, V (j Vi f, Y ( I_+2 ) 

B f Y{ J+2))**3 ‘ ‘ ' ’ _ 

G2=G2+G*XMASS( JJ)«Y{ j + 2}/'RjYtl}» Y( J ), Y(2» ,Y"( j + 1) t YC3) t Y( J+2) j **3 _ 

55 JJ= JJ + 1 J7 

5 CONTINUE " ' 77 ' ' 77. _ "7 

7 F ( I +2 } = FJJ +2 } +G 1— G2 " ~ ~ ' 

: 2 11=1 1+1 _ ~ 

c 7717 _“7 * ’ .7 ~ 77 ' 

C ' ' "7 ‘ ”7 FORCE S 0 N" S UN (_XM ASSll ) ) ' 

c _ ' 7 

'IF IOMPU D.EO) 64,6 5 7_ I 

64' NNL=7_7$NN2 = 1 9 " $ J J=_3 

7G0' TO 66”""'" 7 " 

65 NN 1=4 '?NN2=31__ $_J j£_2 7_ 1 I 

66 "g i=o."" 777 777 7 ' 7 

DC 6 J= NN 1 , MN2, 3 " “ 

"G 1=G 1+G* XNA S'S { J J )"*T yTi ) ~Y ( J ) ) / 1 R ( Y ( ll , Y ( J ) , Y(2 ) , Y i J +T) , Y~( 3"), 

B Y ( J + 2 j ) **3 ’ 71 7_ 

6 JJ=JJ+1 ■ 

f { n=-Gi 



c 


SUNS ACCELERATIONS SET EQUL TO ZERO FOR TEST 
000456 ' Ft I } = 0 « 

" 0004 5 7 IF (CMPL< i) .EQ.l) 6 7,68 "_*" 

'000464 '67 NN 1= 7 $NN2=19 $ JJ=3 " ' 

_000467 ‘ GO TO 69 "H JL ~ .HU ‘ ' 1 ' 

000470 68 NN 1=4 tNN2=3i" $ JJ=2 

_000473 69 G1 = 0. 

000474 ~ DO 7 J=NN1,NN2,3 

1000476 ' G1=G1+G*XMASS( JJ)*(Y{2)-YI J + ll )/R{ Y (1J , Y( J 3 7 Y( 2) Y{ J+l) , Y( 31 * 

VYt J+2) )**3 

' 000521 ' ’7 ' JJ = JJ + 1 

0005 25 F { 2 ) --G i 

0005 27" ~ F { 2 ) = 0 . 

000530 IF (QMPL(I).EO.i) 70,71 J ‘ “ 

000534 70 NN 1= 7 $NN2=19 " iJJ = 3 _ _H _ 

'000537' "GO TO 72 ' • 

^0 005 40^ 7 INN 1= 4 $NN2=3i $JJ = 2 HI 

0005 43 12 G 1=0 . ' ' 

" 000544 ' DO 8 J=NN1»NN2,3 

' 000546 _1 G1=G1+G*XMASS(JJ)*(Y<3)-Yt J+2))/RlY( I),Y( J) »Y12J » Y (J-*- 1 3 ,Y(3) , 

' C Y I J+2 ) ) **3 

Q0~05 71 ~~ 8 JJaJJ-H ' 1H 

"0005 75 F ( 3 ) =— G 1 “ “ 

000577 Fl 3 Y= 0". " 

C ROTATIONAL "MOTION OF"£ARTH 

C" 

C X (34 j = BEt A PRIME'O ~ JH__ 

C X ( 35 ) = BETA PR I ME i™_ _ _ 1 

[ 1 C X(36)=BETA PRIME 2 ' "" " _ 

C X (37) = 8ET A PRIME 3" 

1 c . ' " 

C COMPUTE SETA PRIME AND BET^ PRIME DOT MATRICES ‘ B ^OnVERTEDV 

" 000600" "T=VJDEP-2400000,5 " ■ 

'000602 T= T— 332 82 « + TM 

000605' P 1=3 . 141 592653 58979 " ----- 

000606 TW0PI=2.*Pf 

0006 10' AL=AO+AD*T~ 

'000613 F f 3 7 ) =0 . "" " 

000614 OM { 4 )=0. ■ ' ' 

000615" "DO 109 1 = 1,3 

000616" 0M(i)=0. 

000617 DEL ( I ) = 0 . 

000620" Ft 1 + 33 ) = 0. 



000622 

109 

CONTINUE 








000623 


IF! NOROI 1 ) <, EQ. 0 ) 

GO TO 

1000 






000624 


DO 100 1*1,4 








000626 


BD { I » I ) =V ( 34 ) 








000631 

100 

B ( I , I )=X ( 34 ) 








000635 


B( 1 , 2 ) =X ( 35) 





* 



000637 


B ( 1 , 3 ) -X ( 36) 








000640 


B ( 1 ,4 ) =X( 37 ) 








0C0642 


B { 2 ? 3 )=X ( 37) 








000 643 


B(2 ,4) =-X (36) 








000645 


BI3»4)=X( 35) 







i 

000646 


BD( 1,2>-=V (35) 








000650 


BD ( 1,3) =V< 36) 








000651 


BD ( 1 , 4 ) =V ( 37 ) 








000653 


BD ( 2 t 3)~V (37) 








000654 


BD ( 2 ,4)=~V ( 363 




• 




000656 


8D(3,4)=V(35) 








000657 


00 101 1=2,4 








000661 


JJ^=I~1 








000663 


DO 102 J*1,JJ 






. _ ^ 

_5d„ 

000664 


B ( I , J )--B ( J » I ) 






- 66. 


000670 


B C ( I,J)= -BD{ J , I ) 







k 

000674 

102 

CONTINUE 







h 

000676 

101 

CONTINUE 






gX 

irl 


C 

CALCULATE 0, OMEGA 

1,2,3 

V $CTOR 

,0M 



E 

£ 

000700 


DO 103 1*1,4 






£ 

S 

000701 


DO 104 J* 1 » 4 






§ 

w 

000702 

104 

0M(I)=2.*B( I , J)*V 

(J+33) 

*0M (I ) 




3 

& 

000713 

103 

CONTINUE 









C 

CALCULATE ANGULAR 

R AT--ES 

OF REFERENCE 

AXIS, CAP 

Y 



000715 


FT ( 1 )=0 «. 








000716 


FT ( 2)=X(35)*X(37)- 

-X ( 34 J 

^X(36) 






000723 


FT(3)=X(36)*X(37) 

*X( 34) 

*X { 35 ) 






000727 


FT { 4 ) =X ( 34 ) **2— X { 

35)**2- 

-X(36}**2+X< 37)**2 



• 

000735 


FT ( 2 ) = FT ( 2 ) *2 . *AD 








000740 


FT { 3 ) = FT { 3 ) *2-. *AD 








000741 


FT ( 4 ) =FT ( 4) *AD 








000742 


OM ( 2 ) =0M{ 2 ) +FT ( 2 ) 








000744 


GM(3)=0M(3)+FT(3) 








. 000746 


OM ( 4 )=0M ( 4 ) +Ft ( 4 ) 









c 

CALCULATE MOMENTS 

ACTING ON EARTH PROJECTED Cl\ 

1 BODY AXES Y 




c 

CALCULATE DIRECTION COS 

INES OF 

MOON 

WRT EARTH 

CENTERED AXES, 




c 

LITTLE Y » DEL 










_ 000750 DAL ( 1 ) = X ( 1 3 3 — X { 10) 

000753 DAL ( 2)=X( 14)-X< 11) 

" 000755' DAL (3> = X ( 15 )-X( 12 ) 

C ' ’CALCULATE 8ETA DOUBLE PRIME 

000760 AL= AMOD( AL » TWOP I J 

_ 000 763 C A= COS { A L/ 2 » ) 1. _\L 1 ’ _ _ 

~ 000767 S A= SINUL/2.) 

1000773" " '' B1^CA*X(34)-SA*X(37> 

00100L ' 'B2-CA*X(35)-SA*X< 36) 

001004 B3= S A*X { 3 5 ) +CA#X ( 36 ) 

001010 B4=SA*X(34)+CA*XI 37) 

C CALCULATE ELEMENTS OF ROT AT ION" M AT R IX , C (B ET A DOUBLE PRIME) 

001013 'C(l, 1 ) — B 1 =4= B 1+ B2#B 2—B 3* B3-B 4*B4 

001020 C(l,2)=2.' *( B2*B3+8l*B4) ” 

001024 ' C ( 1 T 3 ) = 2 « *( B2*B4-81*B3) ' 

001030 C ( 2 » 1 ) = 2« * ( B2#B3-B l *B4 ) 

001034 ' C( 2,2 )=BI*B1-B2*B2+B3*B3-B4*84 " 

001040 C( 2,3 )=2.*{B3*B4+B1*82) 

" 001044' ' " C<3,L)=2.*(32*B4+B1*83) "" " 

"‘001050 “ C(3,2)=2.*(B3*B4-B1*62) ‘ 

~001054 C(3,3) = B1*8L-B2*B2-B3*B3 + B4*B4 " 

_0C 1060 ' DO 105 1 = 1,3 ' " 

001062 ' “'DO 106 J=l,3 

001063 106 DELI I >=C ( I , J ) *DAL( J)“+0EL<T) 

0010 74 105 "’CONTINUE 

001076" RRR = SORT { OEL ( i ) **2+DEL (2 j **2+0 EL ( 3 )***2) 

001104 DEL ( 1 )= DEL ( L ) /RRR 

001105 'DEL (2)=DEL< 25/RRR “ 

" 00 110 6’ ' DEL { 3 )=DEL { 3 )/ RRR ’ __ 

C EMIG,EM2G,EM3G, ARE GRAVITY GRADIENT TERMS , 

C FM1,EM2,EM3 ARE ALL OTHER TORQUES " ’ ~ 

001107 _ FM 1 -0 . $ EM2 =0. $ EM3^=0. 

001112’ EM1G=3 » *DEL12)*DFL(3)*XMASS( 5)*ALPE/IRRR**3) 

001120 EM2G=-3." *DEL U ) *DEL ( 3 ) *X MASS ( 5 ) *BE£ /<RRR**3)“ 

001126' EM3G=3. *DEL < l ) *DE L ( 2 )*XMASS { 5 ) *GAF/ < RRR **3 f 

C CALCULATE VALUES’ OF 0,OMEGA1,2,3 DOT VECTOR 

001135 OMD ( L )-0 » " • ” 

_ 0 01136 CM'D ( 2 )•= EM 1G+EM1/E 1 1-ALPE *QM( 3 ) *0M ( 4 ) 

001143 ' OMD ( 3 )=EM2G + EM2/E I2 + BEE*QM< 2 1 *0M (4) 

001151 0MD<4)=EM3G+EM3/EI3~GAE*DM{2 j*0M<3) 

001157 FD(1J=0« 

0 01 r 60 FD ( 2 f=X ( 35 ) * V ( 37 ) Wi'35 ) *X ('37 ) - X“( 3 4 )*V~< 36 ) - V ( 34") * X ( 3 6 ) 

001176 F D ( 2 ) =2 / ’*Ab*FD('2) 



_001201 FD< 3 ) = X ( 36}*V(37KV(36)*X(37)+X( 34) *V(35)+V( 34)*X(35) 

001214 F0 { 3) =2. *AD#FD( 3 ) 

001217 _____ FD(4)=X(34)*V( 34)~X( 35)*V( 35)-X< 36 ) * V ( 36 ) +X ( 37 ) * V { 37 ) 

” 001232' 1 FD ( 4 ) =2 . *AD*FD(4) 

l00 12 35 ~ D(J 107 1=1,4 ' 

001236 DC 108 J = l,4 

_00 1237 11 ' ' >< I+33)=BD( J, I)*(OM( J)-FT( J) ) + B ( J» I ) * (CMD ( J ) -FD ( J ) ) +F ( 1+33 ) 

_ 001267 _ 108 CONTINUE 

J301272 Ft I + 33)=FtI+33)*0.5 _ 

_ 00 12 74 1" 107 CONTINUE _ __ 

”0012 76' ~ "iQCO~ CO NTI NUE " ' • 

C _ ROTATIONAL MOTION OF MOON _ 

C 1 ~ ' , __ 

C X ( 38 )=BETA TRIPLE PRIME 0 

C __ X { 39 5 = BET A TRIPLE PRIME L _ __ 

C XI 40) =BETA TRIPLE PRIME 2 

C XI4D-BETA TRJPLE PRIME 3 


C CALCULATE BETA TRIPLE PR T ME AND BETA TRIPLE PRIME DOT MATRICES 

C BB AND BBD (INVERTED) 

001276 [ I F I NORO 12 ) • FQ • 0 ) GO TO 2000 ' __ 

001277 H'OALI i)=X{ 13 ) — X ( 10 > $ DAL { 2 ) = X (14 )-X { 1 1 } $ DAL ( 3 ) = X ( 15 )-X( 12 )' 


C 

NORMALIZATION OF EULER PARAMETERS 


001307 

XNGRM=X(38)*X(38)+X(39)*X( 39)+X(40) *X< 40} *X( 4 1)*X( 41) 


001316 

XNORM=SQRT< XNORM) 


001320 

X ( 38 )=X( 3 8 ) /XNORM $ XI39) =X( 39) /XNORM $ X(4C)=X< 40) /XNORM 


001327 

X(41)=X(41) /XNORM 


001330 

DO 200 1=1,4 


001331 

BBD( I v 1 ) = VI 3 8 > 

O P, 

001334 

200 BB t I , I ) = X ( 3 8 > 

S?rg 

001340 

BB ( 1 ,2 ) =X ( 39 ) 


001342 

BB ( 1 , 3 ) = X ( 4 0 ) 

Sf' 

001343 

B B < 1 , 4 ) = X ( 41 ) 

*Hri 

001345 

BB ( 2, 3)=X{ 41 ) 



Q 0 1 346 

BB ( 2,4)=-X(40) 


001350 

BBi 3,4)=X(39) 

a § 

001351 

BBD ( 1,2)=V( 39) 

3 sa 

001353 

BBD(i,3)=V(40) 


001354 

BBD { l,4) = V{41 ) 


001356 

8RD(2,3)=V(41) 


001357 

BBO ( 2,4)=— V ( 40 ) 




001361 BBD(3,4)=V(39) 

001362" ' DO 201 1=2,4 

0 0 13 64 _'"JJ=I-1 

‘001366 " ' DO 202 J=1,JJ 

001367" BBII ,J)=~8B< J»I) 

001373 "BBD ( I,J)=~-BBD<j", ff 

0013 77 202~ CONTI NUE 


0014 01 201 CONTINUE 

C “ "‘CALCULATE 0,OMEGAl»2 ,3 VECTOR FOR . MOON 

001403 ~ “ DO 2011 1 = 1,4 

0 014 04 2 01 1' CMM (11=0. 

001407 _ “'DO 203 1 = 1,4' . “ ' 

00 14 10 ^ DO 204 J= 1 , 4 ‘ 

001411 ‘ _ 204 0 M M I I ) = 2. *BB(I,J) *VCJ+37) +OMM(I> 

001422 __ 203 CONTINUE 

001424"' RRR= SORT I DAL 1 1 } **2 +DAL I 2 J **2 +DALI 3 ) #*2 ) 

c ‘ ’ ‘ _ 

C ANGLE "PHI “DEFINED FROM' -PI/2 ‘to" +PI/2 

C 


001433 

CC=DAL ( 1 J /RRR 


00 1434 

CS= D AL { 2) /RRR 


001436 

SPH=D AL ( 3 ) /RRR 


001437 

C PH= SORT (DAL ( 1»**2+DAL( 2) **2)/RRR 


001445 

CL= CC/C PH 


001446 

SL-CS/CPH 


001450 

DAL V I 1 ) =V ( 13 )~V I 10) 


001455 

DAL V(2)=V(14)~V( 11) 


001457 

DAL V{ 3)=V{ 15)-V( 12) 


C R 1 = 0 , R2 = RD0T , R3=RL AMDQTCPH , R4=RPHID0T 

001462 

R 1= 0. 


001463 

R2=CC* DALVtl) +CS*DALV ( 2 ) -i-SPH*DALV< 3) 


001471 

R3= — SL^DAL V ( 1) + CL *DALV<2> ' 


001475 

R4 = -CL*SPH*DALV<1 ) -SL *SPH*DALV { 2 )+CPH* DAL VI 3 ) 


C RT ( 1 ) =0»RT( 2>=~ AMDOTSINPHI, 

C P> T ( 3 ) = PHIDOT , RT I 4)= AMDOTCOSPHI 

001503 

R T { 1 ) =R 1 


001504 

PTI 2)=-R3*$PH/(RRR*CPH) 


001510 

RT ( 3 ) = R4/RRR 


001512 

P T < 4 ) =R3 /RRR 


C CALCULATE AUGMENTED ROTATION MATRIX CIBETA TR«PR IME) ,D 

001513 

DO 220 1=1,4 


001514 

DO 205 J = 1 , 4 


001515 

205 D ( I , J) = 0. 




00 L 5 22 • 220 CONTINUE: 

00 L 5 24 Di 2 ,2 J=X(38 ) **2 + X < 39 1 **2-X ( 40 ) ** 2-X ( 4 il **2 

“001532 _1 '"D(3,3)=X{38)**2-X(39 )**2 + X(40>**2-X(4l)**2 

00L540 D(4,4) = X{38)**2-X{39)**2-X<40)**2 + X{41.3**2 

1.001546] "p<2, 3) = 2. *( X{39 )*X<40 ) *X (38) *X (41 } ) 

0015 53 __ D ( 2 ,4) = 2. * ( X ( 39 ) *X { 41 >-X ( 3 0 ) *X < 40 ) ) 

001560 D( 3 * 2 )=2 • * { X { 39 J*X ( 40 )-X < 38 ) *X ( 4 i ) ) 

'001565 “ D ( 3 > 4 )-2 . * ( X ( 40 ) *X i 41 ) +X {38 ) *X (39 ) ) , 

001572 “ ' 0(4,21=2. * { X ( 39 ) *X (4 1 ) +X (38 1 *X { 40 ) ) 

0015 77 D ( 4 ,3 1 = 2 . * ( X(40 1 *X (41 )-X (38) *X (39 )) __ 

_0 0 1 6 0 4 "_J)0 206 1 = 1,4 _ 

001605 ___ “ _D0 207 J = l,4 “ ‘ ” . “ 

"001606 207 0MM( I )=D( I, J)*RT{ J)+0MM< I) 

"001617“ 206 CONTINUE 

C " CALCULATE MOMENTS ACTING CN MOON 

001621 RRR3=RRR**3 ’ ' 

J301622 PRR2=RRR*#2 „ _ _1L _ 

001623 AMM1=0. $AMM2=0^ $AMM3=0. ' ' ” 

“001626 AMM1G~3 « *0 ( 4 , 2 ) *D ( 3 » 2 ) *XMASS ( 4 1 *ALPM/RR R3 _ 

001633 __ AMM2G=-3. *D { 4 , 2 ) *0 ( 2 ,2 1 * XM AS S( 4) *8EM/RRR3 

_00 1 640 AMM3G=3 . * D ( 3 , 2 ) *D ( 2, 2 1 *X M AS S ( 4 J *GAM /RRR 3 

001645 FAC T- ( 13. 1763965268*3. 1415 926535 8978/180. ) **2 

00165 0 , FACT=FACT*.99Q5*( .0025637252**3) /XMASS14) 

001654 AMM IG=AMM 1G*FACT $ AMM2G-AMM2G*FAC T $ AMM3G* AMM3G*FACT 

C CALCULATE VALUES OF 0»0MEGAi,2,3 DOT VECTOR ' _1 

001657 OMD M ( 1 ) =0 . 

'001660 “ OMDM (2) = AMM1/AMI1+AMM1G-ALPM*OMM(3)*OMM(4> ' ‘ 

"001666 J OMDN(3) = AMM2/ AMI 2+AMM2G+BEM*0MM( 2 ) *0MM ( 4) 

0016 73" OMDM ( 4 )=AMM3/AM I 3+AMM3G-GAM*0MM( 2 1 *QMM{ 3 ) _ 

C CALCULATE VALUES FOR DDT OF RT { I > IE RD( U _1_ 

"2_“ ' C ' CONVERT RT" TO 0 , -LOOTS I NPH I , + PHI DOT , LDOTC OS PH I _ „ 1 


C CALCULATE TIME DERIVATIVE GF ROTATION MATRIX 

001 701 R D { 1 ) =0 o ' J 

001702 RD { 2 )= DALVU) * < R2*SL /RRR2-CL*RT ( 4 >/ ( CPH*RRR > 1 1 

00171 2 RD ( 2) -RD l 2} -DALV (2 )*(R2*CL/RRR2*SL*RT{4i / ( CPH*RR R ) ) 

GO 1 7 22_ ]_RD { 2 )=RD ( 2 ) -( F( 13J- F ( 10 ) )*SL/RRR 

001726 R D { 2 ) = RD(2) + { F ( 14 ) - F ( 11 ) 1 *CL/ RRR ' ~ 

001733 R D ( 2 }— P D ( 2 ) *SPH/C PH +R T ( 3 )*R T (4) / ( CPH*CPH) 

001741 “ ' ' RD(2)=-R0(2) ' ” " 

_001742 RD { 3 3 = DAL V ( 1 )" * {— R T I 2 ) *SL/RRR -J-P2*CL*SPH/RRP2 -rY(3 >*CC/RRRJ 

0017 5 5 "RD < 3 )' '=DA LV< 23* ( .RT ( 2 ) *_ CL./RRR__-RT (3)*CS/RRR " +R2 * SL*SP H/ RRR 21 Y 

_ J 1 R D ( 3 ) ” ~ " “ _ 

001766 RD ( 3 }= DALV (3 ) * (— RT ( 3 ) *SPH/RRR-R2*CPH/RRR2J +RD( 31 



001775 R D { 3 > — ~{ FI 13 )-F { 10 ) ) * CL*SPH/RRR + RDI3) 

'002003' ' RDI3>= -(FI 145-FI im*SL *SPH/RRR+RD< 3 ) 

“002011 ____ RD(3 i= ' ( F{ 15 )-F( 12 i) *CPH/RRR + PDI3) " _ 

' 002016 RD(4)=DAL\/(1J*(RT(2)* CL/IRRR*SPH) ‘ ' 

~ " JL + R 2*SL /RRR2 } " “ 

002025 RD I 4 ) = DALV I 2 ) *" I RT ( 2 ) *SL/ {RRR*SPH ) -R2*CL /RRR2) +RDI4) 

002035 _RD(4)= — (F{ 13)— F( 10) )#SL/RRR + RD(4) 

"002042 " ~RDI4)=(F( 14)- F( 11) )*CL/RRR +RDI4) 

"002047 "" DO 208 1=1,4 ' , 

'002051 " “ DO 209 J=l,4 

"002052 ' ' ~ 2 09 DD( l ,J)=0. 

002057' 208 CONTINUE 

0020 61 " DD ( 2 , 2> — 2 . * ( X I 38 ) *V I 38 ) + X ( 39 )*V ( 39) -X< 40 ) *V { 40 ) -XI 4 1 ) *V( 41 ) ) 

0 020 76 CD { 3 » 3 ) =2 o *{ X( 38 ) *V( 38 ) -X I 39 i *V ( 39) +XI 40 )* V (40 ) -X (4 1 ) * V ( 41) ) 

002113 ^ DD ( 4 ,4) =2 « *1 X 138 ) *V (38 )-X ( 39 )*V ( 39 )-X( 40 ) *V ( 40 ) + XI 4 1 ) *V( 41) ) 

'002130 "DD(2,3)=2. "" * I X(39)*V< 40) +V(39)*XI40) +X (38 i*V (411 + V <38)*XI41 )) 

"002145 " ’ 00(2,48 = 2. * ( X 1 39 ) *V ( 4 L » + V< 39 l*X( 4 1 3 — X ( 38 ) * V { 40 ) -V{ 38 >*X ( 40 ) ) 

00 2162 DD( 3, 2) = 2. * ( X ( 39) *V I 40 ) + V I 39 >*X ( 40 J -X ( 3 8 ) * V ( 41 )- V ( 38 ) *X ( 41 ) ) 

"0021 77 ^ DDI3,4)=2. '*( XI 40 ) *V I 41 ) + V ( 40 )*XI 41) *X { 38 ) *V I 39) + VI 38)*XI 39) > 

002214 DD( 4,2J = 2o " * ( XI 39 > * V I 41 ) + V ( 39 ) *XI 4 1 ) +X (38 }* V 140 ) + V (38 )*X( 40 ) ) 

"00 2231 DDI 4, 3) = 2. ' * I X i 40 > *V 1 41 )"+ V ( 40 j *X< 4.i } -XI 3 8) * V (39 >- VI 38 )*Y< 395 ") 

C 'CALCULATE ACCELERATIONS" 

”002246 ~ DO 2101 i = 'UJtZ 1 II 1_J ~ 

002250 ' FI I+37) = 0. " ' I 

002252' J2101 Ff( U = 0. ' " " 

0022 54 1 " DO 210 1 = 1,4 ' "_L~" _ I" 

002256' J " DO 211 J = l,4 ” ” " I " I_ 

'002257' 211 FT ( 13- ' -D< I , J } *RD( j j-00 (T » j ) *RT ( J ) + FT( I V 

'002300 FTU)=FT( I_) + OMDM ( I ) ‘I 

002303 2ld~ CON T I NUE 


DO 212 1 = 1,4 ____ 

DO 213 J=l, 4 ' "I ‘ 

F I I +37) =F I I +37 5 + BBI J, I)"*Ff (j) 

CONTINUE 

DO 2141 1=1,4" 

FT { I ) = 6. 

DO 2L4 1=1,4 

DO 215 J= 1 , 4 " _ 

FT { I ) =FT { I) + BB(J_7j_)#VXJ + 373* . J__ 

'continue 

'do' 216 1 = 1,4 


~o~o~ 


002344 DO 217 J=l,4 

002345 FI I + 37 ) = F I 1 + 37 )+BBD I J , I> *FT( j 


002-3 53 2 17_C ON T I NUE 

002355 216 CONTINUE 

00235 7"' ' 200Q~R£TURN ^ 

002360 ' "END 
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PROGRAM ANEAMO { INPUT ,OUTP UT , PUNCH ,TA Pfc 2= I NPUT , TAPE3=GUTPUT) 

PROGRAM ANEAMO 

ANALYTIC EARTH AND MOON ROTATION AND TRANSLATION 

THIS PROGRAM CALLS SUBROUTINES LUTH AND LURO PROVIDING RESULTS 
FROM BROWNES LUNAR THEORY AND EC KHARO T* S ROTATIONAL THEORY OF THE 
MOUNTS MOTION AT JULIAN DATES VJIN TU VJF IN INCREMENTS OF VJINC 

THIS PROGRAM CALLS SUBROUTINE EARRQ WHICH USES STANDARD PRECESSION 
-NUTATION FORMULAE AND SIDERAL TIME FORMULAE TO PROVIDE EARTH 
ORI ENT AT ION 

THIS PROGRAM ALSO PROVIDES I/O FUNCTIONS AND CALLS THE 
TRANSFORMATION SUBROUTINES ICOND AND AXANG 

AXANG- CONVERTS A ROTATION MATRIX INTO AXIS AND 
ANGLE OF ROTATION AND T HEN^ INTO EULER 
PARAMETERS 

ICOND- CONVERTS EULER PARAMETERS BETA DOUBLE PRIME 
INTO PARAMETERS BETA PRIME 

ALL ANGLES ARE IN DEGREES, PARALLAX IS IN DEGREES* ALL 
COORDINATES ARE IN KILOMETERS | 

REFERENCES TO MEAN EQUINOX AND ECUPTIC'OF DATE 
REFERENCES 

1) IMPROVED LUNAR EPHEMER IS ( ILEJ 

2) SAO STANDARD EARTH VOl. 1 L966 

_ 1 3) ECKHART A.J. VOL. 70 NO. 7 P.466 

4) WILLIAMS ET.AL. LUNAR PHYSICAL L I BRAT IONS AND LASER RANGING 

MULTI-CASE OPTION " “ I CODE=I " READ NEW TI TLE~ AND' DATES 

0 STOP 


CARD OUTPUT OPTION FOR PARAMETER ESTIMATION 
I P T= 0 NO CARD OUTPUT 

I PT= 1 OUTPUT PHYSICAL L I BRAT I ONS ON CARDS AND OUTPUT INITIAL 

VALUES ON PRINTER __ ' ___ _ 

CARD OUTPUT CONSISTS OF-- ' J ™ 

VJD RHO SIGMA TAU ” ~ 

IN FORMAT (1X.4E19.12) 



000003 


' 000003 


_000003 

000003 

000003 

' 000003' 
000005 
J300007 
OOOQiO" 
"000011 

c 

c 

c 

c 

c 

c 

c 

000017 

000020 

000033 

’ 000046 

000061 
000073 
‘000105 
‘ 000107 

C 

c 

. c 
c 
c 
c 


DIMENSION XEC ( 3 ) » XEQ( 3 ) » VV ( 6 ) *R(3»3),C(3) , BETA ( 4 ) , TT ( 36 ,2) , 

1 CC 131.2) » 

1RR( 31,2) ,NL{ 29) , NLP (29) ,NF (29) ,ND( 29) ,RS(3,3) ,BETDP(4) . SBETA (4) 
1 ,S(3,3),P(3,3),XM0(3,3),BETM1(4),BETP1(4) , BTPDt 4) ,N( 3,3) 

1 ,LKL)K( 14) * LKUA { 14 ) * LKOB ( 14 ) » PKUK ( 1 1 ) , PKOA( 11 ) , PKQB 111) 

2, PL ARG( 23 ) , P LRA TE (23) 

DIMENSION XI ND (21) ,YDEP (21.4) ,TDER(1 ) , FDER ( 1,4) , DER 1( 1,4), 
1DER21 1, 4) ,WK(420) 

2,LATK(26),LATA(26), L AT B( 26 ) , ADA ( 12 ) , A DK ( 12 ) , ADB ( 12 ) 

REAL LKGK, LKDA, LKOB 

COMMON/ LUC UN/ NL, NLP, NF,ND, PL ARG, PIRATE 

COMMON /PERT/ LK.GK, LKOA, LKOB, PKOK, PKQA, PK.OB _ 

l.LATK, LATA, LATB.ADK. ADA. ADB 
P 1=3.14159265358979 
DTK =P I / 180. 

RTD = 180./PI 

NSK IP=0 

MNPTS=2 1 $ NDER=21 $ NCVS=4 $ MMAX=1 $ MDER=1 $ J.W,=-j‘ 

INPUT _ _ 

_ INPUT CONSTANTS FOR ECKHARDTiS THEORY 

J IS THE BETA-GAMMA INDEX- SEE SUBRIUTINE LURO 

J =2 

READ (2,119) (TT( I, J) ,1=1,36) 

READ (2,119) (CC(L, J) »L=1,31) _ _ 

READ (2,119) (RR( M, J ) *M=1 , 31) 

READ (2,119) (PLARG (I ), 1=1,23) 

READ (2,120) (PLRATEt I) , 1 = 1,23) 

DO 6 LL=1 , 29 

6 R EAD( 2.121) N L ( LL) , NL P( LL ) , NF ( LL ) , ND ( LL ) 

INPUT OF SMALL PERIODIC TERMS IN LONGITUDEVLAT ITU'd'E - AND 
PARALLAX FROM THE ILE 

LONGITUDE LIST ( I ALPHA) 

PARALAX LIST (IGAMMA) 

LATITUDE LIST ( IBETA) 



C " INPUT JULIAN DATES. MULTI -CASE" CODE AND OUTPUT' "OPTION 
C 



000124 
' 000126 
000137 
000141 
000143 
' 000154 
' 000156 
000160 


000173 

000175 

000210 

000214 

000232 

000233 


000235 

000245 

000247 

000255 

000260 

000264 

000270 

000276 


000316 

000320 

000322 

000324 

000326 

000331 

000334 


30 


31 


32 


C 

C 

C 

c 


33 

60 


C 

C 

C 


2 


DO 30 ‘1=1.14 

READ ( 2.161 )LKOK( I) * LKOA { I ) » LKOB ( I ) 

CONTINUE 

00 31 1=1,11 

READ (2,161) PKOK(I) , PKOA(I) ,PKOB( I) 

CUNTINUE 

DO 32 1=1,26 

READ (2,166) LATK{ I ) , LATA ( I ) , LATB(I) 


INPUT OF ADDITIVE TERMS IN LONGITUDE ( 1=1,7), NO DE( 1=8, 9) 
AND GAMMAt 1=10, 12) 


DO 33 1=1,12 

READ (2,167) ADK( I ) , ADA { I) . ADB( I ) 

READ (2,900) 

READ (2,100) VJIN,V JINCtVJF, ICQDE, IPT 

IFL AG=0 _ 

VJD=V JIN 


LUNAR THEORY 

CALL LUTH ( V J D , VV , D A, OB , OC , XEC , X EQ ,TU ) ‘ 

I F ( NSKI P.EO.l ) GO TO 20 

IF ( IFLAGoEU. i * AND , IPT.EQ.l) GO TO 20 
WRITE (3,901) 

WRITE (3,900) 

WRITE (3,902) 

WRI TE( 3 ,101) VJD 

WRITE(3»102)VV( 1 ) ♦ VV ( 2) ,VV(3)_,VV(4) ,VV{5),VVjj6) 


C CALCULATE DELAWAY ARGUMENTS ” 

C XL=MOON*S MEAN ANOMALY " _ _ 

C XCAPL=SUN*S MEAN LONGITUDE 

C XLPR=$UN*S MEAN ANOMALY 

C F = MEAN ANOMALY OF MOQN+LUNAR ARGUMENT QF'PERIG EE 

C " ' J. 

20 X L= VV ( 1 ) -VV ( 3 ) 

XCAPL=VV( 1 )-VV( 5) 

XLPR=XCAPL— VV (2) 

F=VV( 1 ) -VV ( 4 ) “ 

XL=AMGD(XL,360. ) [ " 

XLPR=AM0D(XLPR,360. ) 

XCAPL=AMOD( XCAPL, 360. ) 



000337 
000342 
000344 
000347 
0003 53 
000357 
0003 61 
000367 
000372 
000406 
000412 
000424 
000430 
000442 
000446 


000460 

000464 

000506 
000510 
' 000512 
000526 
000527 
000530 
000532 
000536 
' 000544 
' 0005 50 
000564 
000567 
000573 
0005 75 
‘ 000576 
000614 
000616 
000622 
000634 
_000637 
000641 
000645 


F=A MOD ( F * 3 60. ) 

IFIXL.LE.O.) XL = XL + 36 0. 

IF < XCAPL.LE.O. )XCAPL=XCAPL-*-360. 

I F { XLPP. .LE.O. )XLPR=XLPR*360. 

I F ( F.LE.O. } F=F+360 . 

I Fi NSKIP.EQ.l ) GQ TO 21 
IF (IFLAG.EO.l.AND.IPT.EQ.l) GO TO 21 
WRI TE (3,107) 

WRITE( 3.103)XL, XCAPL»XLPR»F 
WRITE(3 ,103) 

WRI TE(3«104)0A»0B»0C 

WRITE (3, 109) 

WRI TE(3»10 5)XEC(1 ),XEC(2)» XEC ( 3 i 
WRITE (3,110) 

WRI TEt 3,106 )XEQ{ 1 ), XE 0(2) ,XEG(3) 

C 

C ECKHARDT* S THEORY FOR LUNAR ROTATION^. __ __ _ ” 

C 

21 CALL EARRO(VJD, R, THETA, S,N,P) ' 

CALL LURO( XL, XLPR.F , VV, P , J ,TAU , Cl , RHO , CC, RR, TT , RS , XMO , XEO*. SLONG 
1 ,S LAT , TU ) 

IF (NSKIP.EO. 1) GQ TO 22 

IF. ( IPT-1 ) 48,49,48 " ■ 

"49 PUNCH 50 » VJD , RHO»C I » T AU ‘ 

IF ( IFLAG.EQ. 0 ) GO TO 48 
GO TO 47 

48 IF (NSKIPoEO. 1) GQ TO 22 
WRITE (3, 113) 

WRITE (3.101) VJD _ ' __ 1 

' WRI TEt 3, 111) "1 “ “ ' 

WR I TE ( 3 * 1 1 2 ) T AU » C I » RH 0 . J “ “ ~ "" ' ' 

CALL AX AN G { R S , DEL, C , BETA) 

WRI TEt 3 ,124) 

DO 7 1 = 1,3 ~ ' ^ 

L=1 

WRITE (3,122) RS ( I , L ) » RS ( I » L+l ) » RS ( I ♦ L+2) 

7 CONTINUE 
WRITE (3, 125) 

WRITE(3,123)( BETAU ), 1 = 1,4) 

22 CALL AXANG ( X MO, DEL , C , BET A ) 

IF(NSKIP.EQ.l) GO TO 23 

WRIT £{3,1 60) “ 

WRI TEt 3, 123) ( BETA ( I) *1=1*4) 



000657 

000661 

000664 


000667 

000670 

000671 

000672 

000674 

000676 

000701 

000702 

000704 

000706 

000716 

000717 

000721 

000724 

000724 

000742 

000744 

000752 

000753 

000754 

000760 


000772 

001000 

7 001007 

001012 

001022 

001026 


DO 10 1=1,4 
10 SBE TA{ I ) =8 ETA ( I ) 

SSLON=SLQNG $ SSLAT =S LAT 
C 

C LOOP TO CALCULATE EULER PARAMETER RATES BETA TRIPLE PRIME BY. 

C NUMERICAL DIFFERENTIATION 

C 

NSKI P=1 
J KJ = 1 
I W=-l 

DE L T= » 0005 „ _ 

TDER ( 1 ) =V JD 6 T STO= TDER ( 1)' 

VJD=VJD-10*DE LT 

GU TO 2 . 

23 XIND( JKJ)=VJD 
DO 24 L = l»4 

24 YDEP( JKJ, L) = 8£TA!L). .. 

JK J= JK J+L 
VJD=VJD+DELT 

IF ( JKJ .GT .21 ) G 0 TQ 25 
GO TO 2 

25 CALL SPLDER I MNPTS , NDER , NCVS , MMAX, MDER, X IND, YDEP , TDER, FDER ♦ 
1DER1.DER2* IW, WK,1ERR) _ __ _ 

DO 27 L=l*4 
27 B TPD t L ) =DER1 ! 1 » L) 

NSK I P=0 
V JD = TST 0 
WRITE! 3,141 ) 

WRI TE ( 3 ,142 ) ( BTPD( I ) » 1 = 1 , 4J 


C EULER PARAMETER TESTS 

T£ST1=SBETA{ 1 )**2+SBETA { 2 ) **2+SBETA( 3 )**2+SBETA <4 ) **2 
TEST2=S BET A! 1 )*BTPD(1 ) +S BETA! 2 5 *BT PD! 2) +S8ETA ! 3 ) *BTPD ( 3 ) + 

1 S8ETA( 4)*STPD(4-) _ : ' 

' WRI TE (3,164) ' ~ 

WRITE (3,165) TE ST1 , TEST2 
WRITE (3,162) 

WRITE (3,163) SSLQN* SSLAT 
C 

C PRECESSION-NUTATION CALCULAtlONS FOR EARTH ORIENTATION 


001036 


40 CALL EARRQIVJD, R, THETA ,S,N»P) 



001042 



001044 



001050 



001056 



001060 



" 001065 



001067 



001072 



001074 



001075 



001113 


12 

"001115 



001121 



001123 



001124 



001142 


13 

001144 



001150 



001152 



001153 



001171 


14 

001173 



‘ 001177 



001201 



001202 



001220 


5 

001222 


41 

"001225 



001227 



001233 



001245 


42 

001251 



001253 


11 

001256 



" 001260 



001264 




C 



C 

LQl 


C 



C 


001276 



001277 



' 001300 




IF (NSKiP.EG.l ) GO TO 41 
WRI TE (3,115) 

WRITE 13,101) VJD 
THETA= THETA*RTD , 

WRI TE ( 3 , 140) THETA 
THETA=THETA*DTR 
WR IT El 3,12 8) 

DO 12 L = l,3 
K=l 

WRITE (3,129) P { L , K ) , P(L*K+1) » P ( L » K+2 ) 

CONTINUE 

WR I TE{ 3,150) 

DO 13 L=1 , 3 
K=1 

WRITE (3,129) ML»K)»N(L,K + 1) »N(L,K+2) 

CONTINUE 

WRI TE( 3, 131 ) 

DO 14 L=1 » 3 
K= 1 

WRITE (3,129) S(l»K)*S(L*K+l) * S ( L» K+2 ) 

CONTINUE 

WRITE (3,124) 

DO 5 L= 1 » 3 _ _____ 

K = 1 

WRITE (3,116) R(L,K) ,R(L,K+1) ,R( L,K+2 ) 


IF (NSKIP.EO.l) GO TO 42 
WRI TE( 3,117) 

WRITE (3,118) I BETA( I), 1 = 1,4“)“ 


THETA, S,N, P) 


DO 11 1=1,4 


IF (NSKIP.EO.l) GO TO 43 
WR I TE ( 3 , 1 3 0 ) 

WRITE (3, 118) (BETA( I) ,1=1,4) 


NUMER-ICAU DIFFERENTIATION 


NSK I P- 1 
JKJ = 1 ' 

DEL T=. 00005 



001 3 02 
001303 
001305 
001310 
001311 
001313 
001315 
001325 
001326 
001330 
001333 
001333 

001351 

001353 

001361 

001362 

001363 

001367 


001401 

001407 

001416 

001421 

001431 

001434 

001436 

001437 

001437 

001441 

001441 

001441 

001441 


001441 

001441 

001441 

001441 


I W=-l _ 

TDER ( 1 ) =VJ D $ TST0=TDER ( 1 ) 

VJD=VJD-10*DELT 
GO TO 40 

43 XI ND ( JK J ) =V JD . 

DO 45 L = 1 * 4 

45 YDEP( JKJ*L)=BETA( L) 

JK J- J KJ+ 1 

V JD=VJD+DELT 

IF { JKJ.GT.21 ) GO TO 44 

GO TO 40 _ . 

44 CALL SPLDERt MNPTS * NDER* NCVS* MMAX *MDER» X IND* YDEP » TDER* FDER » 

1 DER1.DER2* IhUWK* IERR) 

DO 46 L=1 ♦ 4 

46 BETDP(L)=DERltl*L) . . . 

NSK I P=0 

” VJD=T ST 0 „ 

WP.ITE<3*127> 

WRITE13 »126) (BETDP (I), 1 = 1.4) 

C EULER PARAMETER TESTS 

TEST1=SBET A( 1 )**2+SBETA{2)**2+SBETAl3)**2*SBETA(4)**2 
TE ST2=S3ETA( 1 ) * BETD P( 1 > +SBETA (2 ) *BETDP( 2) +SBETA (3 ) *8ETDP{ 3) 

1 +SBETA(4)*BETDP(4) 

WRITE (3*164) 

WRITE (3*165) TE ST1 » TEST2 _ 

47 IF (VJD.GE.VJF) GO TO 3 _ 

VJD=V JD+V J I NC _ ^ ' _ 

IFLAG=1 _ 

GO TO 2 

3 IF ( ICODE.EO.l ) GO TO 60 
50 FORMAT C1X,4E19.12) 

100 FORMAT (3E20. 10*215) 

101 FORMAT {* *»44X»*JULIAN DATE= *.E19.12*/ ) 

102 FORMAT {* *»*MUON= * » F10.5 *5X ♦ #GAMMA= **F10.5*5X* * GAMMA PRIME = *, 
1F10 . 5* / ** *,*OMEGA= # ♦ F 1 0. 5 » 5X* * D= * , F 10 . 5 , 5X * *0B L I QUI TY= #,F10.5* 
2///) 

107 FORMAT (50X ,* DELAUNAY ARGUMENTS* * / ) 

103 FORMAT (* *,*L= *, F 10 .5 * 5X , *CAP L= * , F10.5* 5X**L PRIME= *»F10.5*5X 

lt*F= *»F10. 5, ///) _ _ J L 

108 FORMAT { SOX * * ECL I PT IC LONG. AND LAT.*»/) 

104 FORMAT (* * ,*LONGITUDE= *, FlO. 5*5X* *L AT ITUDE=*»F10.5»5X,* PARALLAX 



001441 

001441 

001441 

001441 

001441 

001441 

001441 
001441 
001441 
"001441 
001441 
001441 
001441 
001441 
001441 
001441 
001441 
001441 
001441 
001441 
001441 
001441 ' 
001441 
001441 
001441 
001441 
001441 
001441 " 
001441 
001441 
001441 
001441 
001441 
001441“ 
•001441 
001441 
001441 

001441 • 
001441" 
001441 


1=*.F10. 5,///> 

109 FORMAT (50X ,* RECTANGULAR ECLIPTIC COORDINATES*./} 

105 FORMAT {* *.#X= *,F 12.4, 5X. *Y= * , F 12 . 4, 5X , *Z = *»F12.4.///) 

110 FORMAT (50X ,* RECTANGULAR EQUATORIAL COORDINATES*. / ) 

111 FORMAT (50X,*PHYSICAL L I BRAT I ONS*. / ) 

106 FORMAT (* *»*X= * » F 1 2 «4»5X»*Y= * , F 12 . 4. 5X . *Z= *,F12.4,///} 

112 FORMAT (* * » * LONG „= * , F 12 . 9 . 5 X , *NODE= * , F 12. 9 . 5X , * INC L I NAT ION= *, 
1F12.9,5X,*BET A GAMMA INDEX = *,I3,///> 

113 FORMAT (*1*.50X»*MOONS ORIENTATION*.///) 

119 FORMAT ( 7F10 . 5 ) 

120 FORMAT I5F15.5) 

121 FORMAT (412) 

117 FORMAT (50X.*EULER P ARA METERS -ME QEQ TO BODY*,/) 

124 FORMAT (30X,*D IRECT ION COS INES { MEQEQ5 0 TO BODY)*,/) 

122 FORMAT (20X , 3 ( E19 .12 , 5X ) ) 

125 FORMAT ( / / , 50 X , *E UL ER PARAMETE RS-MEQEG TO BODY*,/) 

123 FORMAT ( 20X » 4 ( E 19 .1 2 , 5X ) , / ) _ ‘1 ‘ 2__' 

128 FURMAT I SOX , * PRECES S I ON MATRIX*,/)’ 

115 FORMAT C *1* , 5 OX , *EARTH ORIENTATION*,/) 

116 FORMAT (20X,3{E19»12,5X) ) 

118 FORMAT {20X,4(E19.12,5X)»/ ) 

127 FORMAT ( // , 50 X ,*EUL ER PARAMETER RATES-REF TO BODY*,/) 

126 FURMAT { 20X.4 ( El 9. 1 2 . 5X) , //) " 

129 FORMAT ( 20X , 3 { El 9 .1 2 , 5X ) ) 

130 FORMAT (50X,*EULER PARAMET ERS -RE F . TO BODY*,/) 

131 FORMAT ( / *50X ,*SPIN MATRIX*,/) 

140 FORMAT (44X,*S IDEREAL T IME=*, E 19 . 12 , / > 

141 FORMAT ( SOX , * EULER PARAMETER RATES REF. TO BODY*,/) 

142 FORMAT { 20X , 4 ( El 9.1 2 , 5X ) , / // ) 

150 FORMAT ( / , SOX ,*NUTATI ON MATRIX*,/) 

160 FORMAT 150X.*EULER PARAMETERS REF. TO BODY*,/) 

161 FURMAT (3E20.12 ) 

162 FORMAT (50X,*EARTHS SELENOGRAPH IC COORDINATES*,/) 

163 FURMAT ( 20X , * LONG=* , E 19 . 12 , 10X, *LAT=* . E 19.J.2, / ) 

164 FURMAT ( / , 50X ,*EULE R PARAMETER TESTS*, / )~ 

165 FORMAT { 40X , E 22. 15 * 5X , E 22 . 15 ) ' _ ' ' __ 

167 FURMAT ( 2E 10. 3, E20. 12 ) ““ ” ” " 

166 FORMAT ( E10.3 .E10.3 . E20 .12 ) “ “ " ” _ " " 

900 FORMAT (8QH 

1 ) 

901 FORMAT (*1*) 

902 FORMAT (* *,50X,*LUNAR ORBIT*,// ) 

STOP 





SUBROUTINE LU RO{ ARGLv ARGLP, ARGF , VV , P» J ,T AU ,S IG« RHO.C, R, T, RM, XMG, 
1 XSOtSLONGtSLAT *TU) 


C 

C THIS SUBROUTINE PROVIOES THE PHYSICAL L I BRAT ION IN LONGITUDE (TAU), 

C ' “'IN NODE { IS IGMA) » AND IN INCLINATION (RHO) FROM ECKHARDT+S LUNAR 

C L IBRATI ON T AB LE S-RE F : TH E MOON VI P264. 

C ' ADDITIVE AND PLANETARY TERMS ARE INCLUDED PER REF-' 

C LUNAR PHYSICAL L I BRAT IONS AND LASER RANGING, W ILL I AMS , ET .A L > 

C 

C_ INPUT 

C T ( I , J ) =TAU COEFFICIENTS 
C I =TERM NO. J=COEFPICI ENTS { BETA) 

C ' J = 1 BET A=0* 0006268 GAMMA=0 .0002300 (GT 1 ARC-SEC) 

C J=2 BET A=0. 00063 GAMMA=0.0002 +ADDI TI VE/PLANETARY 

I__I_. C ' _ TERMS ( G T .3. ARC-SEC) 

II C ~~ ARGL=MEAN ANOMALY OF THE MOON= ' ' ” 

C' ARGL P=MEAN ANOMALY OF THE SUN= 

C ARGF=MOQN-QME GA “ 

C AP.GD=MEAN ELONGATION OF MOON FROM SUN(D) 

C C ( I » J )-SIGI COEFFICIENTS 

C R ( I , J ) = RH 0 COEFFICIENTS 

C 

000025 ” DIMENSION T l 36 , 2) ,C ( 3 i, 2 i » R( 3 1, 2 ) , NL ( 29 ) , NLP ( 29 ) , NF( 2 9) . ND (29 ) 

000025 " DIMENSION RM( 3 , 3) , EC { 3, 3) , RR { 3, 3) , VV ( 6 ) , PLARGI 23 ) ,PLRATE(23) 

000025 , DIMENSION X MO ( 3 , 3 ) , XE 0 ( 3 ) , XMU 3 , 3) , P ( 3, 3 ) , XE05 { 3 ) ' 

000025 "" COMMON/ LUC ON/ NL , NLP , NF» ND, PLARG , PLRAT E 

000025 ' ' ARGD= VV { 5 ) " ' ' . 

000026 “ IF ( J.EO.l )XI=5521. 5 

000031 ■” IF ( J.EQ.2 >XI=5549.3 I I __ „ _I~“ _ 

000034' P I = 3. 14159265358979 " " 

000036 ' T WOP 1 = 2 .*PI 

000037 ' RTD=1 80 ./PI " 

000041 DTK=P I / 180. 

000042 TTT =TU* 36525 « 

000044 ARGL = AR GL’J'DT R 

000045 ARGL P = ARGL P*DTR ' 

000046 ARCF=ARGF*DTR 

000047 ARGD=ARGD*DTR 

000050 " TAU = 0. $ SIG I =0 « $ RHG=0. 

000053 DO 1 1 = 1,13 ~ " “ 

000055 ARG=NL( I ) *ARGL+NLPII) *ARGL P+NF { I )*ARGF+ND( I)*ARGD 



000071 
QUO 106 
000107 
0001 U 

000131 

000145 

000147 

000150 

000 170 
000205 
000206 
000212 
000213 
000216 
000220 
000234 
000250 
000265 
000273 
000277 
000301 
000306 
000317 
000322 
000325 
000331 
000334 
000341 
000345 


000351 

000356 

000365 

000367 

000375 

000402 

000406 

'000410 

000411 


1 TAU=TAU+T{ It J )* SIN1ARG) 

I D I V=13 

DO 2 1=1,3 

APG=NL { I+IDIV )*ARGL+NLP { I+IDIV) *ARGLP+NF( I + ID IV )*ARGF+ND{ I+I DIV) * 
AARGD 

2 SIGI=SIGI+CtI , J)* SIN(ARG) 

I DI V= ID IV + 8 

DO 3 1=1,8 

AP.G=NL ( I + IDIV )*ARGL+NLPt I + I DIV J*ARGLP+NF{ I+IDIVJ^ARGF+NDt I+IDIV)* 
BARGD 

3 RHG=RHU+R{ I, J)* COS(ARG) 

DO 50 1=1,23 

AA=PLAP, G( I ) +PLRATEI I ) *TTT 
AA= AA*T WOP I 
AA= AMOD ( AA , TW OP I ) 

IF (AA.LT.O.) AA= AA + T WOP I 

TAU=TAU+T< 1+1 3, J)*S I N t A A ) ' ' * 

SIGI=SIGI+C( I+8»J)*SIN£AA) 

50 RH(J = RHO +R £ I +8 » J ) *COS £ AA ) 

TAU= TAU/3600. $ RHQ=RH0/36Q0. $ SI G= I S I G I /X I ) *RTD 
TH=RHO+XI/3600. 

PSI=VV£ 4)+SIG 
PHI=180-+VV( 1 )-PSI+TAU 

PSI ,=AMODt PSI ,360.) $ P H I = AMOD { P H 1 , 36 6 « J $ t H= AMO D £ TH , 360 . } 

IF (PSI.LT.O. ) PSI=PSI+360. 

IF <PHI.LT.0.i PH 1= PH 1+360 . 

IF (TH.LT.O.) TH=TH+360. 

PSI=PSI*DTR $ PHI=PH1*DTR $ TH= TH*DTR' ' ' ' " 

C PH= COS < PM I ) $ SPH= SIN(PHI) _ 1 

CPS= COS(PSI) $ SPS= SIN(PSI) 

CTH = COS(TH) S STH= SIN£ TH) 

C 

C RM- ECLIPTIC OF DATE TO BODY ROTATION MATRIX 

C 

RM{ 1,1) =CPH*CPS-SPH*CTH*SPS ' ' " 

RMt 1,2)=-CPH*SPS-$PH*CTH*CPS $ R M £ 1,2) = -RM ( 1 ,2 )~ 

RMt 1»3)=-SPH*STH 

RMt2,l)=CPS*SPH+CPH*CTH*SPS $ RM( 2, 1 ) =-RM{ 2,1 1 

RM£2»2)=-SPS*SPH+CPH*CTH*CPS 

RM£ 2,3) =C PH*STH $RM (2 ,3) =-RM( 2,3) 

RMt 3, 1) =-SPS*STH 

RM t 3, 2 1 = C PS* STH 

RMt 3» 3) =CTH 



000413 
000421 
000424 
000435 
000465 
000467 
000470 
000475 
000477 
000500 
000501 
000502 
0005 17 
000521 


000523 

000524 

000525 

000530 

000532 

000550 

000552 

000554 

000555 

000556 

000560 

000573 

000575 

000576 


000602 

000612 

000614 

000615 

000617 

000625 


VV( 6)=VV(6)*DTR 

EC{ 1, 1 ) =1 . $ ECl 1, 2 )=0. $ EC { 1» 3 J = 0. 

EC ( 2» 1 ) =0 . $ EC ( 3 1 1 ) =0 . $ EC ( 2t 2 ) = C0S(VV{6)) 

EC { 2 1 3 ) = SlNlVVioJ) $ EC{3,2)=- SIN{VV<6)) $ EC(3»3)= C0S(VV(6)) 
DU 12 I =1 1 3 
DO 13 L = 1 » 3 
13 RR{ I, 0 = 0, 

12 CONTINUE 
DO 10 1=1,3 
DU 11 L= 1 , 3 
DO 15 K=l,3 

15 RRtI,L)=RM{I,K)*EC(K,L)+RRCI»U 
11 CONTINUE 
10 CONTINUE 
C 

C MULTIPLY MATRIX RR BY P TO REFER ANGLES TO MEAN EQUATOR AND 
C_ EQUINOX OF 1950 o 0 __ __ 

C ' ' RM- MEQEQ50 TO BODY ROTAT ION' M ATRIX' 

C 

DO 30 1=1,3 
DO 31 L=l»3 
RM( I, L) = 0. 

DO 32 K=l,3 

32 RM( I,LI=RR(I*K)*P(K,L')+RM{I,L') 

31 CONTINUE 
30 CONTINUE 

DO 33 1=1,3 
XEQ5 ( I ) =0. 

DO 34 K=l,3 

34 XEQ5U)=P(K,n*XEQ{K)+XEQ5U') 

33 CONTINUE 
DO 35 1=1,3 

35 XEQII) = XEQ5(I ) 

C 

C CALCULATION OF EULER PARAMETERS FOR ROTATION FROM JJPCASE__2_ 

C' " TO LOWCASE Z FOR USE AS INITIAL CONDI T ION IN' RIGEM~_ 

C " 

RRR= SORT (XEQ ll»**2+XEQ(2) **2+XEQ ( 3 ) **2 ) 

CC= XEQ( 1) /RRR 
CS= XEQI2 ) /RRR 
S PH = XEQ{ 3 } / RRR 

' CPH= S0RT(XEGm**2+XEQ(2)**2T/RRR 

SL= XEQI2 ) /{RRR*CPH) 



000630 


C L = XEQ(l)/ (RRR*CPH) 

000631 


XM1( 1*1 )“-CL*CPH 

$ 

XM 1 1 1 , 2 ) =SL $ XMlt 1,3)=~CL*SPH 

000636 


XM1 { 2*1)=-CPH*SL 

$ 

XM1{ 2»2)=-CL $ XM1 ( 2, 3)=— SL^SPH 

000641 


XM 1 { 3 * 1 ) =-S PH 

$ 

XM 1(3*2) =0. $ XM1 (3 , 3 )=CPH 

000644 


DO 20 1=1,3 



000651 


DO 21 K=i * 3 



000652 

21 

XMO ( I , K )=0 • 



000660 

20 

CONTINUE 



000662 


DO 22 1=1,3 



000663 


DO 23 K=1 * 3 



000664 


DO 24 L = 1 , 3 



000665 

24 

XHO { I , K ) =XHO { I * K.) +R M ( I, L ) *XM1 ( L , K) 

000704 

23 

CONTINUE 



000705 

22 

CONTINUE 




000710 

000716 

000724 

000732 

000735 

000736 


C 

C CALCULATION OF EARTHS S ELENOGRAPH I C COORDINATES 
C __ 

SS= SORT CXMOt i.IJ**2+XM0(2fl)**2> 
SLUNG=ATAN2 (XMQ (2, 1 ) , XM0( 1,1)1 
SLAT=ATAN2 ( X MG ( 3, 1 ) , SS ) 

SLONG=S LONGER TD $ SLAT =SLAT*RT 0 

RETURN 

END 



c 

c 

_____ c 

c 

c 

c 

c 

c 

c 

_ ” c 

c 

c 

__ c 
c 

“ 1__ c 

c 

c 

c 

c 

c 

1 "lie 

c 

000013 


000013 

000013 

000013 

000027 

000030 

000033 

000034 

000035 

000040 

000042 

000043 

000046 


000105 


SUBROUTINE LUTH ( V J D, VV t OA , GB , OC , XEC » XEO t TU ) 

THIS SUBROUTINE PROVIDES AN APPROXIMATE VERSION OF BROWN* S LUNAR 
THEORY {REF. I L E ) 

INPUT 

THE CALLING PROGRAM SHOULD PROVIDE THE JULIAN DATE (VJD) 

VARIABLES 

VV1 -MEAN LONGITUDE OF MOON* MEASURED IN ECLIPTIC FROM MEAN EQUINOX 
OF DATE TO MEAN ASC. NODE OF LUNAR ORBIT THEN ALONG ORBIT (DEGREES) 
{ MOON) 

VV2-SUN *S MEAN LONGITUDE OF PERIGEE (DEG) (GAMMA) 


VV3 = MEAN LONGITUDE OF LUNAR P ER I GEE ,MEAS . IN ECLIPTIC FROM MGAU EQUINOX 
OF DATE TO MEAN ASCENDING NODE OF LUNAR ORBIT THEN ALONG ORBIT 
(DEG) (GAMMA PRIME) 

VV4 =LQNG. OF MEAN ASC. NODE OF LUNAR ORBIT ON ECLIPTIC MEAS. FROM 

MEAN EQUINOX OF DATE (DEG) IOMEGA) 

VV5=MEAN ELONG. OF MOON FROM SUN XEC ( 1 , 2, 3 J = ECL IPT IC RECTANGULAR 
COORDINATES (DEG) (D) 

0A» G8»0C= MOUNTS ECLIPTIC LONGITUDE LATITUDE PARALLAX ( DEG . ,DEG . , DEG. ) 


XEO ( 1,2,3)=EQUAT0RI AL RECTANGULAR COORDINATES 

V6=0BLI QUITY OF THE ECLIPTIC (DEG) (OBLIQUITY.) _ * L.1_Z 

DIMENSION XEC (3 ) , XE 0 { 3 ) , VV { 6 ) ,L KOK { 14 ) , LKOA( 14)»LK0B{ 14), ‘ 

1 PKOK ( 11) , PKOA{ 11 ), PKOB ( 11 ) “ " ’ 

2 , LATK(26) ,LATA(26) ,LATB(26), ADAU2J. ADBt 12) ,ADK( 12) " 

REAL LKOKt LKGA « LKOB 

COMMON/PERT/ LKOK»LKOA»LKQB» PKOK »_PK0A, PKOB ' 

1, LATK, LATA, LATBvADK »ADA, ADS 

DMSTR(D«VM»S)= 3. 141592653589793 ' /180. ”*( D+( VM+S/60." )/60._ ) 

PI=3. 141592653589793 

TU= (VJD -2 415020. ) Z36525. 

TU2=TU*TU 



TU3=TU2*TU _ ___ 

C0R1 = 13 36 . *2. *Pl““ 

COR 2= 11. *2 . *PI 

COR3=5. *2. *PI 

C 0R4= 1236. *2. *P I 

1 Vl= DMSTRC270. ,26. *3.69 ) + ( COR 1 + DMSTR(307. *52. *59.31 

_10})*TU - DMSTR ( 0. *0. ,4.08' )*TU2 + DMST_R(0 0. ,.0068 _)* 

~ 2TU3 ' ' 

V2= DMSTRI 28L . ,13. ,15. ) + DMSTR (0 • ,0. ”,6139.03 )#TU ' 



000 1 44 
000204 

0002 44 ' 

0003 05 

000343 
000346 
000347 
000351 
□00353 
0003 54 
0003 56 
000357 
000362 
000365 
000370 
000373 
000376 
00040L 
000404 
000407 
000413 
000417 
000423 
_000434 __ 
000441 

C 

C 

C 

000447 

'000456 

000460 

000467 

000477 

00050l__ 

000510 


14- DMSTFUO . ,0. ,1.63 )* TU 2 + DMSTR ( 0. ,0. ,.012 )*TU3 

73= DMSTR ( 334. ,19. ,46.75 ) + (C GR2 + DMSTR U 09 . ,2. ,2.52 >) 

1 *TU~ DMSTR{0. ,0. ,37.17 } *TU2-DMSTR { 0. ,0. ,.045 ) *TU3 
V 4= DMSTR (259. ,10. , 59.79 ) - (CDR3+DMSTR( 134 . ,8. ,31.23 ) 

1)*TU + DMSTR 1 0. ,0. ,7.48 )*TU2+ DMSTR 1 0 . , 0. ,.008 >*TU3 

75= DMSTR{ 350 . ,44. ,15.65 ) + (CGR4+ DMSTR(307. ,6. ,51.18 

1))*TU- DMSTR ( 0. ,0. ,5.17 )*TU2+ DMSTR { 0. ,0. ,.0068 )*TU3 

3 76= DMSTR ( <23 . ,27. ,8.26 } - DMSTR ( 0 . ,0. ,46.845 )*TU~ 

1 DMSTR ( 0 • ,0. ,.0059 >*TU2 + DMSTR! 0. ,0. ,.00181 )*TU3 

V1=V1*180. /PI 

72= V2*l 80 . /PI ' ' 

V3=V3*1 80 . /PI 
V4= 74*1 80. /PI 
75=7 5*180. /PI 
76=76* 180. /PI 
R360=360. 

7 1=AM0D (71,360. ) ' " ] 

V2=AM0D l V 2 , 360. ) ' ' " ' 

73=AM0D (V3 ,360. ) 

74=AM0D (74,360. ) 

75=AM0D(V5,36O. ) 

76 = A MOD <76,360. ) 

I Ft VI *LE. 0* ) VI = 71 +3 60 . __ 

I F { V2.LE.0, ) V2=V2+360. 

I F t V3.LE.0.) V3=V3+360. 

IFIV4.LE.0.) 74=74+360. 

IF(V5.LE.O.) V5=V5+360, 

VV{1)=V1 S7V( 2)=V2 $VV ( 3 }=V3 $ VV(4)=V4 $ VV<5)=V5 $VV{6)=V6 

V 1= V 1*P I / 1 80. $ V2=V2*P I /I 80. $ V3=V3 *P I / 180 . 

V4=V4*P 1/180. $ V5=V5*PI/180.'$"V6=V6*PI/180V 

CALCULATION OF ADDITIVE TERMS ’ ' _..1~ ” 

TE1=. 53733431 - ( 10 1049 82. E-l 2 ) *TU*36525 . " 

1 +191»6-16*TU*TU *36525. **2 ‘ 

TEL=TE1*2. *PI • ~ 1 

ATL = 14. 27 * SIN(TEl)*PI/(3600. *180. )' 

TE2=. 71995354 -{ 1470942 28. E- 12 ) *TU*3 6525 _ 

1 +43.E-16*TU*TU*36525« **2 ' J_ 

TE2=TE2*2. *PI 

AT0=95 . 96* SI N( TE2 ) *P 1/ ( 3600* *180. 0) ~ 

T E3 =. 48398132 - ( 1472 69 147 . E— 12) *TU*36525 _ . 

1 +43.E-16*TU*TU*36525 . **2 



000520 TE3=TE3*2. * P I 

000522 " ATU 1 = 1 5 .5 8 * $ IN ( T E3 ) *P I / { 3600 . *180. 

000531 T E4=. 71 9953 54- { 14709422 8. E-12 ) * TU*36525 . 


1 + 143. E-16)*TU*TU*36 525. **2 


000541 

“000543 

,000552 

’ 000562 
000564 
000573 
'000577 
~ 0006 04 
000612 
000626 
000627 
' 000635 
“000651 
000652 
000660 
' 000674 
000677 


_ TE4=TE4*2 . *P I 

ATL 1=7.261*5 INI T54 ) *P I / ( 36 00. *1 8 0. ) 

TE5=. 52453 688- t 1471 62675. E-12 )*TU*365 25. 
1 +( 43.E-16)*TU*TU*36525.**2 
TE5=TE5*2.*PI 

ATU2 = 1. 86*SINITE5)*PI/13600.*180) 

ATL2=0. $ AT 03=0 . $ DGC=0. 

DO 13 1=1.17 

ADARG= ( ADAt D+ADDi I >*TU*36525.) *2.*PI 

13 ATL2=ATL2+ADK ( I )*SIN{ ADARG} 

DO 14 1=8,9 

ADARG= I ADA ( I) +ADB ( I ) *TU*36525. J *2.*PI 

14 ATO 3=AT03 + S IN ( ADARG )* AD K ( I ) ' 

DO 15 1=10,12 

ADARG= ( ADA ( I ) +ADB ( I ) *TU*36525 . ) *2. *P l" 

15 DGC=DGC +ADK ( I } *COS (ADARG ) 
VI=V1+ATL+ATL1+ATL2 
V4=V4+AT0+AT01+AT02+AT03 


C 

C CALCULATION OF PERIODIC TERMS 

C 1 


000704 

VA = V 1 

~V3 

000706 

V3=V1 

-V5 

000710 

VC=VB- 

-V2 

000712 

VD=V1 

-V4 

000714 

V5S= 

SINIV5) 

000716 

V5C= 

COS( V5 J 

000720 

V 52=2*V5 

000723 

V52S 

= SINIV52) 

000725 

V52C= 

COS ( V52 ) 

000727 

V 54=4*V5 

000732 

V54S= 

SI N ( V5 4 ) 

000734 

V54C= 

C0SIV54) 

000736 

VAS= 

SINI VA) 

000740 

VAC = 

COSIVA) 

000742 

VA2=2*VA 

000745 

VA2S = 

SIN{ VA2) 

000747 

VA2C= 

C0SIVA2) 

000751 

V A3=3*VA 



C00754 

VA 3S = 

S I N{ VA 3 J 

000756 

V A3 C- 

CGS< VA 3) 

000760 

VB S= 

S I N ( VB ) 

000762 

V BC= 

COS( VB ) 

000764 

VCS = 

S IN { VC ) 

“000766 

VCC = 

COS( VC ) 

000770 

vos= 

S INI VD 1 

000772 

VDC= 

COS ( VO ) 

000774 

VU2=2*VD 

000777 

VD2 S= 

SI N ( VD2) 

001001 

V D2 C= 

COS ( VD2 ) 

"001003 

Vlll = 

VA-V52 

001005 

V11S = 

S INC VI 11) 

001007 

V11C= 

COStVlll) 

001011 

VI 21 = 

VA2-V52 

001013 

V 12S= 

S I N ( V 1 2 1 ) 

001015 

V 13 1= 

VA+-VC-V52 

001020 

V13 S= 

SI N ( VI 31 ) 

001022 

V 13 C = 

COS ( VI 31 ) 

001024 

V141 = 

VA-5-V52 

001026 

V14 S= 

S I N ( V 1 41 ) 

001030 

V 14 C= 

C0S(2*VA+V52) 

001036 

V14CC 

= COS ( VI 41 > 

‘001040 

V15 1= 

VC-V52 

001042 

V15S 

= SINIV151) 

001044 

V 1 5C = 

COS( V151) 

001046 

V16 1= 

VA-VC 

001050 

V16S= 

S I N ( V 1 6 1 ) 

001052 

V16C= 

C0S{ Vlfcl) 

001054 

V17 1= 

VA+VC 

001056 

V17S= 

SIN ( VI 71 } 

001060 

V 1 7 C = 

COS I VI 71 ) 

001062 

V 18 1 = 

VD2-V52 

001064 

VI 8 S= 

SINC V181 J 

001066 

V 19 1 = 

VA+VD2 

001070 

V19 S= 

S I INI { VI 91 ) " 

001072 

V 21 1= 

VA-VD2 

001074 

V21C= 

COS { V2 11 ) 

001076 

V21S 

= SINIV211) 

001100 

V3 11 = 

VA-V54 

' 001102 

V31 S= 

SI N ( V311) 

001104 

V31 C= 

COS ( V3 1 1 ) 

‘ 001106 

V41 1 = 

VA2-V5*4 



001112 
001114 
00111 ? 
001121 
001125 
001125 
001127 
001131 
001133 
001135 
001137 
001141 
001144 
001146 
001151 
001153 
001155 
001157 
001161 
001163 
001166 
001170 
001172 
001174 
0011 77 
001201 
001203 
001205 
001207 
001211 
001213 
001215 
001217 
001221 
001224 
001226 
001232 
' 001234 
001237 
001241 
_001243 
001245 
001247 


V41 S= SIN ( V4 1 1 ) 
V5U = VA-VC-V52 
V51S= SUUV511) 
V61 1=VC+V52 
V61S= S I N { V611J 
V71 1 =VA +VD 
V71 S= SIN(V7in 
V81 1 = VD-VA 
V01S= SI N( V8 11 } 
V9 1 1= VD-V5 2 
V 91 S= SI N ( V9 1 1 ) 
M 22 1 = VD +V52-V A 
V22S= SIN(V221) 

V231-VD+VA-V52 
V23S - S I N ( V231 ) 
V24 L=VD+V52 
V24 S= SI N ( V241 ) 
V25 1=VA 2 + VD 
V25 S= SIN( V2 51 ) 
V261=VD-V52~VA 
V26S= S I N ( V 2 6 1 ) 
V271=VD-VA2 
V27S= SI N ( V271 ) 
V28 1=VC +VD-V52 
V28 S= SINIV281) 
VV1 1=-V 16 1 
VV1S= SlN(VVll) 
V291-V611-VA 
V29 S= SINW291) 
V441=V151-VA 
V44S= S I N { V441 } 
V45 1=-V21 1 
V45S= SIN( V451 ) 
V54=V5#4 
V54 S= SI N ( V54) 
All=V5*2+VA2 
A1 1 S= SIN (All) 
A12=VA-VC+V52 
A 12S= S I N ( A 1 2 ) 
A13=VA~V5 
A13S= SINIA131 
A 14=VC+V5 
A14S= SINIA14) 



001251 

A15=VA3-V52 



001253 

A15 S = 

SINIA15) 

. 

001255 

V61C = 

COS ( V6 1 1 ) 


001257 

V 12C = 

COS { V 1 21 ) 


001261 

V4 1 C= 

COS { V4 1 1 ) 

. 

001263 

A12C = 

C0SU12) 


001265 

V51 C = 

COS (V511J 


001267 

AA1 = VA2-VC $ AA2=VA2+VC 

$ AA3=V A3- V52 

001275 

AA 1 C= 

CUS( AA1 ) $ AA2C= 

COS ( AA 2 ) $ AA3C= C0S(AA3J 

001303 

AA4 =VC + V5 $ AA5=VA + V5 $ 

AA6= VD2-V52 

001311 

AA4 C = 

COS ( AA4 ) $ AA 5C = 

C0S(AA5) $ AA6 C= C0S1AA6) 

"001317 

QA= 22639.5 *VAS -4586.463 *V11S +2369.912 *V52S +769.016 * 


1 VA2S -668*146 *VCS -411.608 *VD2S -211.656 *V12S -205.962 

2 * V 1 3 S +191.953 * V14S -165.145 *Vi5S +147.607 *V16S -125. 

3 154 *V5S -109.673 *V17S -55.173 *V18S -45.099 *V19S+39,528 

4 * V21S -38.428 *V3 LS +36.124 *V A3S -30.773 *V41S +28.475.. * 

5 V5 1S -24.420 *V61 S _ „ 

6 + 13.902 * V54 S + 14.387 * A11S + 14.577 * A12S + 18.609 * 

7 A 1 3S + 18.023 * A14S - 13.193 * A15S ‘ 

DO 10 1=1,14 

ALKARMLKOAC U+LKUBi I >*TU*36525. 1*2. *PI 

10 0A=UA+LKGKl I ) * SIN(ALKAR) 

0 A= (OA/ 3600 . 1*PI/180. 

0 A= OA+V 1 

(JO 3422.70 +186.5393 *VAC +34.3117 #V11C +28.2333 *V52C 

1 +10.1657 *VA2C +3 .0861 *V14CC+1.9178 *V15C +1.4437 *V13C 

2 +1.1528 *V16C -.9781 *V5C -.9490 *V17C -.7136 * V21C +.6215, 

3 *VA3C +.6008 *V31C 

4 + .2607 * V54C - .3 * V61C - .3997 * VCC + .2833 * V14C 

"' 5 - .3039 * V12C + .3722 * V41C + .2302 “*"A12C - .2257.' * 

1 V51C .... 


001411 

001412 

001420 

001433 

001436 

001440 


001512 

001527 

001531 

001537 

001552 

001554 

001557 

001560 


0C = 0C+» 1268 *AAlC-ol038 *AA2C-.1187 *AA3C+.1494 *AA4C 

1 1093*AA5C-.1052*AA6C 

DO 11 1=1,11 

PKARG=(PKOAU J+PKOB ( I)*TU*36525. )*2. *PI 

11 OOOC+PKQK(I}* CQS(PKARG) 

0C=0C*(1. -4. 6747E-5 ) 

0C3={GC*0C*0C )/{6. *206265. **2) 

UC=0C+0C3 

0BB=- 112.79 *V5 S+2 373. 36 *V52S+192.72 *V14S+226 09. 07 *VAS 

1 -4578.13 #V11S -38.64 *V31S +767.96 *VA2S -152.53 *V12$ -34. 

2 07*V41 S+50. 64 *VA3S -25.10 * V61 S ' -126. 98 *VCS -165.06 *V15S 

3 -115.18 *V17S -182.36 *V13S -23.59 *V29S -138.76 *VV1S -31.70 



001630 
001631 
00 1641 
_001 6 55 
001660 

'001675 
001676 
001706 
C01722 
001724 
001736 
001746 
001752 
001753 
001756 
'001760 
001765 
001767 
001771 
001772 
002005 
002010 
0020 12 
002013 
002021 
002027 
002035 


002 043 
002045 
Q 02 050 
002052 
0020 56 
002060 


002062' 

002064 


4 #V44S -52.14 *V18S~85.13 #V45S 

DU 12 1=1,22 

LARG= (LATA { I J +LAT B ( I 3 #TU#36525. ) #2.#PI 
12 0BB=0B8 + LATKU)*S!N(LARG) 

0 BB = 0BB#PI/(180. #6 0. #60. ) 

0883= -526.069 #V91S #44.297 #V23S #20.599 #V81S -30.598 #V26S 

1 -24.649 #V27S -22.571 #V28$ 

DO 16 1=24,26 

LARG= (DATA ( I 3 +LATB ( I )#TU*36525. 3 #2.#PI 
16 0BBB=06B8 + LATKU)*S IN(LARG) 

S=VD+OBB 

SS= SI N t S3 $ SS3 = S IN { 3 .#S ) $ SS5 = S INI 5.#S) 

LARG=ILATA(23)+LAT8(233#TU#36525.}#2.#PI 

COK=1.+2.708E~6#139»978#DGC 

G1=18519«75#C0K 

G2=6 . 24 1 #C0K##3 

_ G3=.004#C0K#*5 . . . 

SSA=G1#LATK( 2 33 *S IN ( LARG ) . 

SSB=G2#SSA/G1 

SSC=G3#SSA/G1 

SSD=SSA/G1 _ 

0B= SSA# SS + SSt3#SS3 #SSC#SS5#0B88 

OB=OB#PI/( 180. #60. #60. ) _ 

OC=OC#PI/ ( 180. #60. #60. ) 

2 HP= 6378 . 16 /OC _ . 

0AS= SIN(OA) " ’ ” _ “ ' . “ ’ . 

0 AC = COS t OA ) _ ... > 

OBS= S I N { OB ) 

OBC= COS ( OB 3 _ _ 

c 

C CALCULATION OF ECLIPTIC AND MEAN EQUINOS OF DATE COORDINATES, 

C _ 

H PX=HP#OBC#OAC _ 

HPV=HP#OBC#OAS _ 

_ HPZ=HP# DBS ’ * 1 

XEC ( 1 3 =HPX $ XEC ( 2 )=HPY' $ XEC{3)=HPZ 

V6S= SIN(V6) .... -• • 

V6C= COS( V63 

C ' __ 

C CALCULATION OF MEAN EQUATOR AND EQUINOX OF DATE COORDINATES, 


RADGX=HPX 

RADGY=HPY#V6C-HPZ#V6S 



0020 70 
002072 
002100 
002106 
002107 
002110 
002113 
002116 
002121 
002122 


RADGZ={ HPY)*V6S+HPZ*V6C 

XEO { 1 )=RADGX $XE 0 { 2 ) =RADGY $ XE0(3)=RADGZ 

GA=(JA# 1 80 . /PI 

0B=0B*18U. /Pi 

00=00*180. /PI 

OA = AMOD{QA s R360) 

GB= AMOD (0 B » R3 60 ) 

00= A MOO { OC » R36 0 1 

RETURN 

END 



SUB ROUT INE 6 ARRO { VJ Dt R , T HETA , S t N, P ) 


OOGOll 
OOOOI1 
oooo n 
000012 
000014 
0000 16 


c 

C THIS SUBROUTINE UTILIZES NEWCOMS ANALYTICAL EXPRESSIONS TO CALCULATE 

C . THE PRECESSION MATRIX «P * THE NOTATION MATRIXrN? AND THE SIDERAL 
C ROTATION MATRIXES. INPUT IS THE JULIAN DAT E » V JD . OUTPUT CONSISTS OF 

C THE DIRECTION COSINES RELATING THE X AND W AXIS SYSTEMS OF THE 

C REFERENCE IN THE FORM X=(SNP)W 

C 

C REF: SAO STANOARD EARTH 1966 

C 

REAL N(3»3i 

DIMENSION S ( 3 * 3 ) « P { 3» 3 ) » R ( 3 »3 ) 

PI = 3. 141592 65 358979 
TWO PI* 2. *PI 
DTR=P I / 180 . 

RTD= 1 80 . /PI 


C 

C 

_C 

000017 
’ 000021 
1 0000 23 
000024 
Q00030 
“ 000032 
‘ 000036 
0000 43 
000046 _ 

“0000 51 ' 
000054 
000057 


000 1 24__ 

000127 

000135 

000143 

000150 

C 

C 

C 


CALCULATE MODIFIED JULIAN DATE AND .S 

XMJD=VJD-2400000.5 

T=XMJD-33282.0 

T2=T*T 

ARG 1*1 12.1128 • -.052954 *T)*DTR _ __ 

ARG2=ARG1*2. " _ 

ARG3=( 280.0812 +.985647 *T)*DTR*2. 

AP G4- ( 64. 3824 +13.176398 *T)*DTR*2. 

ARGl= AMUDt ARG1 ,TWUPI > 

ARG2= AMGDI ARG2.TW0PI ) _ 

ARG3= AMUD { ARG3 ♦ TWO P I ) 

ARG4= AMODI ARG4,TWUPI ) “ 

THETA* (100.075542 +360.985647348 *T+.29E-12*T2-4. 392E-3* 

ASIN( ARG1) +.053E-3* SIN( ARG2) 325E-3* SI N C ARG3 ) - . 05 E— 3* 

BS IN ( ARG4))*DTR .... 

THETA=AMOD{ T HET A » TW OP I ) _ _ ’ 

"CTH* COS{ THETA) 

' STH= SI N( THETA) ' " " “ 

SU«1)=CTH $ S(2» 2 )=CTH $ S ( 1 ♦ 2 ) =STH'' i S(2,1)=-STH ' " “ _ 

S ( 1 s 3 ) = 0. $ S(3 » 1) *0. $ S ( 3 9 3 ) *1 . $ S(2*3)*0.. J _S (3t2)*0. 

^CALCULATE. NUTATION MATRIX, N 


_v .. 


000156 


Rl* ARG1 



000157 

000161 

000162 

000205 

000227 

000251 
000255 
000265 
000275 
000306 
000311 
0003 14 
000316 
000322 
000325 


000327 

000332 

000334 

000337 

000400 

000443 

000460 

000523 

000565 

000602 

000620 

000635 


000644 

000646 

000647 

000652 

000653 

000654 

000672 


R2=ARG3 . . ; 

R3= ARG4 

DEL MU = - 7 6.7E-6* SIN(Rl) +.9E-6* SIN(2*R1) -5.7E-6* SIM f R 2 }- 
1 .9E-6* S I N < R3) 

DEL NU= -33.3E-6* SIN(Rl) +.4E-6* SIN(2*R1> -2.5E-6* SI N< R2 ) _-.4E-6 
" 1 * SIN( R3> , 

DELEP= 44.7E-6* COS(Rl)- .4E-6* C0S(2*R1) +2*7£-6* C0S(R2) 

1 +.4E-6* C0S(R3) 

CNU=COS (DELNU ) * SNU=S IN ( DELNU } 

CMU= CDS! -DELMU ) $ SMU= SIN(-DELMU) 

CEP= CCJS { -DEL EP ) $ SEP= SIN(-DELEP) _ _ „ 

N( 1,1)=CNU*CMU 6 N(1.2)=CNU*SMU $ N{1,3) =-SNU 

N<2.1)=CMU*SEP*SNU-SMU*CEP 

N(2 *2) = SMU*SE P*SNU+CEP*CMU 

N ( 2 * 3 ) = SEP#CNU ' ... 

N(3»1)*CMU*SNU*CEP+SEP*SMU 

N ( 3 * 2 ) = SMU# SNU^CEP-SEP^CMU __ 

N(3 *3)=CEP*CNU 
C 

C CALCULATE PRECESSION MATRIX »P 

C . . 

XKAP=0. 063107 *T*DTR/36O0. 

UMEG=0. 063107 *T*DTR/3600. ' 

XNU=0. 0548757 *T *D TR/3600 . 

P(1,I)=- S I N{ XKAP ) * S I N ( OM£G ) -s- COS(XKAPJ* COS(OMEG)* COSUNUJ 
P ( 1 ♦ 2 ) = - COS ( XKAP) ^ SIN(OMEG)- SIN(XKAP)* -COS( OMEG) * COS(XNU) 

P { 1 , 3 ) = - COS ( OMEG ) * SIN(XMU) 

P { 2 » 1 ) = SINfXKAP)* COS (OMEG) + COS(XKAP)* SIN(OMEG)* COS(XNU) 

P ( 2 * 21= COS ( X KAP ) * COS(OMEG)- SIN(XKAP)* S I N( OMEG )*..COS < XJMU) 

P (2 * 3 1 = - S IN ( OMEG ) # SIN(XNU) ’ ..... 

P t 3 * 1 ) = COS ( XKAP) * SIN(XNU) 

P (3 *2) = - SINtXKAPl# SIN(XNU) 

P 1 3 » 3 1 = COS(XNU) 

C 

C CALCULATE FINAL ROTAT ION_ MATRIX *_R 

c ” 

DO 1 K=l*3 ‘ ’ _ _ _ 

DO 2 J=l,3 
R(K»J}=0. 

DO 3 1=1*3 " _ _ 

D04 L = 1 * 3 

R(K*J)=R{K*J) +S ( K* I )#N( I»L)#P(L*J) 

4 CONTINUE 



000675 3 CONTINUE 

000677 " 2 CONTINUE 

000701 I CONTINUE 

000703 RETURN 

000704 END 




SUBROUTINE ICONDt VJD, BETAt 


THETA, S.N.P) 


C 

C THIS SUBROUTINE PROVIDES THE TRANSFORMATION FROM BETA 

C DOUBLE PRIME EULER PARAMETERS TO BETA PRIME PARAMETERS 

C 

C BETA ENTERS AS BETA DOUBLE PRIME 

C BETA RETURNS AS BETA PRIME 

C 


000011 

"000011 
000011 
000013 
000015 
000016 
000020 
0000 22 


000023 
0000 27 
000032 
000034 
000036 
000040 
000044 
000050 
000054 


000060 
000065 
000071 
0000 72 
0U0073 
000104 


DIMENSION BET A{4> ,B (4,4 ) , BT ( 4, 4 ) , BET AS { 4 ) , OM ( 4 ) , B ETADt 4 ) , B ETDP ( 4) 
1 * S ( 3 « 3 } , P t 3 » 3 ) » 01 ( 3 > 

REAL M { 3 , 3 ) 

XMJD=VJD-2400Q00.5 
T=XMJD-33282 . 

P 1= 3. L41592 65358979 
TWO P 1=2 , *PI 

DTR=P I / 180 . ' i_ 

RTD = 180. /PI 
C 

C BETA MATRIX, 6 

C 

ALPHA={100. 075542 +360,985647348 *T)*OTR 

ALP HA=AMODt ALPHA, TWOPIJ 
ALP HA= ALPHA/2 . 

C A2= COS( ALPHA ) 

SA2= S I N( ALPHA ) 

B 1 1 , 1 ) =CA2 $ 8 ( 1 , 2 ) = 0 . $ B ( 1 ♦ 3 ) =0 « $ Btl,4)=~SA2 

B(2,U = 0, $ B ( 2, 2 ) =C A2 $ Bt2,3)=-SA2 $ 8{2,4)=0. 

B ( 3 , l ) = 0. $ B{ 3 , 2 i =S A2 $ B(3»3)=CA2 $'Bt3,4) = 0. 

B (4 » 1 ) =S A2 $ B{ 4, 2 ) = 0 « $ B(4,3)=0 o $ B(4,4J=CA2 
C 

C TRANSPOSE B TO GET BETA INVERSE, BT 

c ' 

DO 5 1=1,4 ' __ 

5 BETAS ( I ) = BETA { I ) __ J ' ______ '' 

00 1 1 = 1,4 ’ 

DO 2 J= 1,4 
2 BT { I, J)=B( J,I ) 

1 CONTINUE 
C 

C CALCULATE BETA PRIME 



000110 

000112 

000113 

000125 

000127 

000130 


BET A ( I ) =0 . 

00 4 J=1 » 4 

4 BET A( I ) =BET Ai I)+BTI I, J) *BETAS ( J } 
3 CONTINUE 
RETURN 
END 



SUBROUTINE AXANGIR, DEL* C. BETA ) 


C 

C THIS SUBROUTINE CALCULATES THE AXIS AND ANGLE OF REVOLUTION FOR 

C_" AND ROTATION MATRIX. R. IT ALSO PROVIDES THE EULER PARAMETERS 

C "BETA 
C 

C REF : KORN AND KORN 

C 


000007 DIMENSION R ( 3 . 3 ) . C ( 3 ) , B ET A (4 J 

000007 PI=3. 14159265358979 

000010 " TP=R(1»1)+R{2.2)+R{3.3) 

000013 ' CDE L= ( TR-1 . )/2. 

000016 SDE L= SORTll* -CDEL*CDEL) 

000022 ' DEL= ATAN2I SDEL, CDE L ) 

C 

__ ■ C DEL IS SUPLLIED IN RANGE ZERO TO PI 

C " 


000031 

000035 

000041 

000045 


000051 

000052 

000056 

000060 

000062 

000067 

000075 

000103 

000111 

000112 


C(1}=(R(3»2)~R(2»3) 1/12. *$DEL ) 

C(2) = (R{i»3)“RI3»l) ) / ( 2 . * S D E L ) 

C(3)=IR{2tl)-Rll»2) )/(2. *SD6U 
C(l)=-C(ll 6 C l 2) --C { 2J $ C ( 3 )=-C( 3 ) 

C ..... 

C‘ SEQUENCE TO KEEP CALCULATED ROTATION AXIS GENERALLY ALIGNED" WITH 
C BODY ROTATION AXIS 

C “ 

IF {C(3).GEc0. } GO TO 1 

C(l)=-CilJ S C I 2 J = -C ( 2) $ C(3)=~C( 3) 

DEL=2 . *PI-DEL ‘ 

1 DEL=0EL/2. 

C CALCULATE EULER PARAMETERS "" 

BETA! 1) = COS( DEL ) 

BETA(2)=C(1)* SIN(DEL) 

BETAI3)=C(2>* SIN(DEL) 

BET A{4J-C{3)# SIN(DEL) __ __ 

RETURN" 

END 



